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A reflection of our ultimate understanding of a complex system is our ability to control 
its behavior. Typically, control has multiple prerequisites: it requires an accurate map 
of the network that governs the interactions between the system’s components, a quan¬ 
titative description of the dynamical laws that govern the temporal behavior of each 
component, and an ability to influence the state and temporal behavior of a selected 
subset of the components. With deep roots in nonlinear dynamics and control theory, 
notions of control and controllability have taken a new life recently in the study of com¬ 
plex networks, inspiring several fundamental questions: What are the control principles 
of complex systems? How do networks organize themselves to balance control with func¬ 
tionality? To address these here we review recent advances on the controllability and 
the control of complex networks, exploring the intricate interplay between a system’s 
structure, captured by its network topology, and the dynamical laws that govern the in¬ 
teractions between the components. We match the pertinent mathematical results with 
empirical findings and applications. We show that uncovering the control principles of 
complex systems can help us explore and ultimately understand the fundamental laws 
that govern their behavior. 
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I. INTRODUCTION 

To understand the mechanisms governing the behavior 
of a complex system, we must be able to measure its state 
variables and to mathematically model the dynamics of 
each of the system’s components. Consequently, the tra¬ 
ditional theory of complex systems has predominantly 
focused on the measurement and the modeling problem. 
Recently, however, questions pertaining to the control of 
complex networks became an important research topic 


in statistical physics (Cornelius et al. 20I3| Liu et al. 

2011a 

Nepusz and Vicsek| 2012[ Ruths and Ruths |2014 

Sun and Motter 

2013 Yan et al. 

2012 

). This interest is 


driven by the challenge to understand the fundamental 
control principles of an arbitrary self-organized system. 
Indeed, there is an increasing realization that the design 
principles of many complex systems are genuinely deter¬ 
mined by the need to control their behavior. For exam¬ 
ple, we cannot divorce the understanding of subcellular 


networks from questions on how the activity or the con¬ 
centrations of genes, proteins, and other biomolecules are 
controlled. Similarly, the structure and the daily activity 
of an organization is deeply determined by governance 
and leadership principles. Finally, to maintain the func¬ 
tionality of large technological systems, like the power 
grid or the Internet, and to adapt their functions to the 
shifting needs of the users, we must solve a host of control 
questions. These and many similar applications have led 
to a burst of research activity, aiming to uncover to what 
degree the topology of a real network behind a complex 
system encodes our ability to control it. 


The current advances in controlling complex systems 
were facilitated by progress in network science, offering 
a quantitative framework to understand the design prin- 


ciples of complex networks (Albert and Barabasi 

2002 

Barabasi and Albert 

1999 

IDorogovtsev et al.. 2008 

]Mii^ 

et al. 2002 Newman 

2006| Toroczkai and Bassler 

2004 

Watts and Strogatz 

1998 

I. On one end, these advances 


have shown that the topologies of most real systems share 
numerous universal characteristics. Equally important 
was the realization that these universal topological fea¬ 
tures are the result of the common dynamical princi¬ 
ples that govern their emergence and growth. At the 
same time we learned that the topology fundamentally 
affects the dynamical processes taking place on these 


networks, from epidemic spreading ( 

Cohen et al. 

.. 

Pastor-Satorras and Vespignani 2001 

I to synchroniza- 

tion ( 

Nishikawa et al. 

2003 

Wang and Slotine 

20051. 


Hence, it is fair to expect that the network topology of a 
system also affects our ability to control it. 


While the term “control” is frequently used in numer¬ 
ous disciplines with rather diverse meanings, here we em¬ 
ploy it in the strict mathematical sense of control theory, 
a highly developed interdisciplinary branch of engineer¬ 
ing and mathematics. Control theory asks how to in¬ 
fluence the behavior of a dynamical system with appro¬ 
priately chosen inputs so that the system’s output fol¬ 
lows a desired trajectory or final state. A key notion 
in control theory is the feedback process: The difference 
between the actual and desired output is applied as feed¬ 
back to the system’s input, forcing the system’s output 
to converge to the desired output. Feedback control has 
deep roots in physics and engineering. For example, the 
centrifugal governor, one of the first practical control de¬ 
vices, has been used to regulate the pressure and distance 
between millstones in windmills since the 17th century 
and was used by James Watt to to maintain the steady 
velocity of a steam engine. The feedback mechanism re¬ 
lies on a system of balls rotating around an axis, with a 
velocity proportional to the engine velocity. When the 
rotational velocity increases, the centrifugal force pushes 
the balls farther from the axis, opening valves to let the 
vapor escape. This lowers the pressure inside the boiler, 
slowing the engine (Fig. [^. The first definitive math¬ 
ematical description of the centrifugal governor used in 
















































































3 



FIG. 1 Feedback control. A centrifugal governor represents a 
practical realization of a feedback process designed to control 
the speed of an engine. It uses velocity-dependent centrifu¬ 
gal force to regulate the release of fuel (or working fluid), 
maintaining a near-constant speed of the engine. It has been 
frequently used in steam engines, regulating the admission of 
steam into the cylinder(s). 


Watt’s steam engine was provided by James Maxwell in 
1868, proposing some of the best known feedback control 
mechanisms in use today (Maxwell 18681. 

The subsequent need to design well controlled engi¬ 
neered systems has resulted in a mathematically sophis¬ 
ticated array of control theoretical tools, which are today 
widely applied in the design of electric circuits, manu¬ 
facturing processes, communication systems, airplanes, 
spacecrafts and robots. Furthermore, since issues of reg¬ 
ulation and control are central to the study of biological 
and biochemical systems, the concepts and tools devel¬ 
oped in control theory have proven useful in the study 


of biological mechanisms and disease treatment (Iglesias 


and Ingalls 2009 Sontag 2004). For example, feedback 


control by transcranial electrical stimulation has been 
used to restore the aberrant brain activity during epilep¬ 
tic seizures (Berenyi et al. 2012). 


Modern control theory heavily relies on the state space 
representation (also known as the “time-domain ap¬ 
proach”), where a control system is described by a set 
of inputs, outputs and state variables connected by a set 
of differential (or difference) equations. The concept of 
state, introduced into control theory by Rudolf Kalman 
in 1960s, is a mathematical entity that mediates between 
the inputs and the outputs of a dynamical system, while 
emphasizing the notions of causality and internal struc¬ 
ture. Any state of a dynamical system can then be rep¬ 
resented as a vector in the state space whose axes are 
the state variables. The concept of the state space was 
inspired by the phase spaee concept used in physics, de¬ 
veloped in the late 19th century by Ludwig Boltzmann, 
Henri Poincare, and Willard Gibbs. 


For a nonlinear dynamical system, we can write the 


state space model as 


x(t) = f(t,x(t),u(t);0) 
y{t) = h(t,x(t),u(t);0) 


(la) 

(lb) 


where the state vector x(t) G represents the internal 


state of the system at time t, the input vector u(t) G 
captures the known input signals, and the output vector 
y(t) G captures the set of experimentally measured 
variables. The functions f(-) and h(-) are generally non¬ 
linear, and 0 collects the system’s parameters. Equations 
(la) and are called the state and output equations, 
respectively, and describe the dynamics of a wide range 
of complex systems. For example, in metabolic networks 
the state vector x(t) represents the concentrations of all 
metabolites in a cell, the inputs u{t) represent regulatory 
signals modulated through enzyme abundance, and the 
outputs y{t) are experimental assays capturing the con¬ 
centrations of a particular set of secreted species or the 
fluxes of a group of reactions of interest. In communica¬ 
tion systems x(f) is the amount of information processed 
by a node and y{t) is the measurable traffic on selected 
links or nodes. 

A significant body of work in control theory focuses on 
linear systems ( Kailath[ 1980), described by 


x(t) = A{t) x(t) 
y{t) = C{t)x{t) 


B(t)u(t) 

D(t)u(t), 


(2a) 

(2b) 


where (2a) and (2b) represent so-called linear time- 
varying (LTV) systems. Here, A(t) G is the 

state or system matrix, telling us which components 
interact with each other and the strength or the na¬ 


ture of those interactions; B(t) G 


IxM 


matrix; C(t) G 


IxN 


is the input 


is the output matrix; D(t) G 


is feedthrough or feedforward matrix. In 
case A(t), B(t), C(t) and D(t) are constant matrices. 


(2a) and (2b) represent a linear time-invariant (LTI) sys¬ 


tem, which is the starting point of most control theo¬ 
retical approaches. Note that since we typically know 
u(t) and D(t), we can simply define a new output vector 
y{t) = y{t) — D(t)u(t) = C{t) x(t), allowing us to ignore 
the D(t)u(t) term. 


Many nonlinear systems like (la lb) can be linearized 


around their equilibrium points, resulting in an LTI sys¬ 
tem. For example, in stick balancing, a prototypical con¬ 


trol problem (Luenberger 1979), our goal is to balance 


(or control) the stick in the upright position using the 
horizontal position of the hand as the control input u(t). 
This mechanical system has a natural state space repre¬ 
sentation derived from Newton’s second law of motion. 
Consider a stick of length L whose mass M is concen¬ 
trated at the top. Denote the angle between the stick 


^ For a more realistic case, treating the stick as a rigid body of 
uniform desnity, see ([Stepan and KoUar[ [2000^. 
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FIG. 2 (Color online) Two mechanical systems whose natu¬ 
ral state space representation with linear time-invarian (LTI) 
dynamics can be derived from Newton’s laws of motion, (a) 
The goal of stick balancing, a simple but much studied con¬ 
trol problem (also known as the inverted perdulum problem), 
is to balance a stick on a palm. Redrawn after ( Luenberger[ 


19791. (b) A rocket being thrust upward. The rocket ascends 


from the surface of the earth with thrust force guaranteed by 
the ejection of mass. Redrawn after (Rugh 19931. 


and the vertical direction with 0{t). The hand and the 
top of the stick have horizontal displacement u(t) and 
x{t), respectively (Fig. [^. The nonlinear equation of 
motion for this system is 

L0{t) = g sin 0{t) — u{t) cos 0{t), (3) 

where g is the gravitational constant and 

x{t) = u{t) + L sin 0{t). (4) 

When the stick is nearly at rest in the upright vertical 
position (0 = 0, which is an equilibrium point), 0 is small, 
hence we can linearize Eqs. (§ and obtaining 

x{t) = jr[x{t) - u{t)]. (5) 

Using the state vector x(t) = (x(t), with velocity 

v(t) = x{t), and assuming y(t) = x{t), we can rewrite the 
state and output equations in the form of an LTI system 



0 1 


' 0 ' 

x(t) = 

i 0. 

x(t) -b 

9 

L_ 


i{t) 


yit) = 


1 0 


c(t). 


(6a) 

(6b) 


This form allows us to perform linear controllability 
analysis. Indeed, as we show in Sec |II.B[ the linearized 
system (6a I is controllable, in line with our experience 


that we can balance a stick on our palm. 

Linearization of a nonlinear system around its normi- 
nal trajectory {x*(t),u*(t)} generally leads to an LTV 


system. Consider the motion of a rocket thrust upward, 
following 


m{t)'h{t) = rh{t)ve — 'm{t)g, (7) 

where m{t) is the mass of the rocket at time t and h{t) 
is its altitute. The thrust force m{t)ve follows Newton’s 
third law of motion, where m{t) denotes the mass flow 
rate and Ve is the assumed-constant exit velocity of the 
exhaust (Fig. If we dehne the state vector x(t) = 
{h{t),v{t),m{t))^ with velocity v{t) = h{t), the control 
input u(t) = rhit), and the output y{t) = h(t), we have 
the state-space representation 


Xi{t) 


X2(t) 

X2{t) 

= 

u(t)v^ 

xs{t) 


. . 


y{t) = xi{t). (9) 

The state equation ([^ is clearly nonlinear. Let’s con¬ 
sider its linearization around a nominal trajectory that 
corresponds to a constant control input u*{t) = uq < 0, 
i.e. a constant mass flow rate. This nominal trajectory 
follows xl{t) = Ve[{mo/u(j + t)ln(l -|- uot/mo)] — gt^/2, 
X2{t) = Ve ln{l + Hot/mo)—gt, x^it) = mo+uot, where mo 
is the initial mass of the rocket. By evaluating the partial 
derivatives and at the nominal tracjectory, 

we obtain the linearized state and output equations in the 
form an LTV system 


/ 

'o 

1 

0 


0 

^s{t) = 

0 

0 

— UQVe 

^sit) + 

Ve 

(mo+uot)^ 

mo+uot 


_o 

0 

0 


1 

II 

1 

0 

0 X5(t), 




us{t) 


( 10 ) 

where the deviation variables xs{t) = x(t) — x*(t), 
us(t) = u{t) - and ys{t) = y{t) - y*{t) = X 5 (t). 

Notwithstanding our ability to design such well- 
controlled systems as a car or an airplane, we continue to 
lack an understanding of the control principles that gov¬ 
ern self-organized complex networked systems. Indeed, 
if given the wiring diagram of a cell, we do not under¬ 
stand the fundamental principles that govern its control, 
nor do we have tools to extract them. Until recently 
the degree of penetration of control theoretical tools in 
the study of complex systems was limited. The reason is 
that to extract the predictive power of (laI and (lb), we 
need (i) the accurate wiring diagram of the system; (ii) 
a description of the nonlinear dynamics that governs the 
interactions between the components; and (iii) a precise 
knowledge of the system parameters. For most complex 
systems we lack some of these prerequisites. For exam¬ 
ple, current estimates indicate that in human cells the 
available protein-protein interaction maps cover less than 
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20% of all potential protein-protein interactions (Sahni 


et al. 2015); in communication systems we may be able 


to build an accurate wiring diagram, but we often lack the 
analytical form of the system dynamics f (x(t), u(t); 0); 
in biochemical reaction systems we have a good under¬ 
standing of the underlying network and dynamics, but 
we lack the precise values of the system parameters, like 
the reaction rate constants. Though progress is made 
on all three fronts, offering increasingly accurate data 
on the network structure, dynamics, and the system pa¬ 
rameters, accessing them all at once is still infeasible for 
most complex systems. Despite these difficulties, in the 
past decade we have seen significant advances pertaining 
to the control of complex systems. These advances in¬ 
dicate that many fundamental control problems can be 
addressed without knowing all the details of equations 
(la) and (lb). Hence, we do not have to wait for the 


description of complex systems to be complete and ac¬ 
curate to address and understand the control principles 
governing their behavior. 

Graph theoretical methods have been successfully ap¬ 
plied to investigate the structural and the qualitative 


properties of dynamical systems since 1960’s (Yamada 


and Foulds, 1990). This raises a question: Can the re¬ 


cent renaissance of interest in controlling networked sys¬ 
tems offer a better understanding of control principles 
than previous graph theoretical methods? To answer 
this we must realize that the current interest in control 
in the area of complex systems is driven by the need 
to understand such large-scale complex networks as the 
Internet, the WWW, wireless communication networks, 
power grids, global transportation systems, genome-scale 
metabolic networks, protein interaction networks and 
gene regulatory networks, to name only a few (Chen 


2014). Until the emergence of network science in the 


21th century we lacked the mathematical tools to char¬ 
acterize the structure of these systems, not even mention¬ 
ing their control principles. The non-trivial topology of 
real-world networks, uncovered and characterized in the 
past two decades, brings an intrinsic layer of complexity 
to most control problems, requiring us to rely on tools 
borrowed from many disciplines to address them. A typ¬ 
ical example is the structural controllability problem of 
complex networks. Structural control theory developed 
in 1970’s offered sufficient and necessary conditions to 
check if any network with LTI dynamics is structurally 


controllable (Lin 1974). Yet, it failed to offer an efh- 


The goal of this article is to review the current ad¬ 
vances in controlling complex systems, be they of biolog¬ 
ical, social, or technological in nature. To achieve this we 
discuss a series of topics that are essential to understand 
the control principles of networks, with emphasis on the 
impact of the network structure on control. The review 
is organized around several fundamental issues: 

(i) Controllability. Before deciding how to control a 
system, we must make sure that it is possible to control 
it. Controllability, a key notion in modern control the¬ 
ory quantifies our ability to steer a dynamical system to 
a desired final state in finite time. We will discuss the 
impact of network topology on our ability to control com¬ 
plex networks, and address some practical issues, like the 
energy or effort required for control. 

(ii) Observability. As a dual concept of controllability, 
observability describes the possibility of inferring the ini¬ 
tial state of a dynamical system by monitoring its time- 
dependent outputs. We will discuss different methods to 
identify the sensor nodes, whose measurements over time 
enable us to infer the initial state of the whole system. 
We also explore a closely related concept — identifia- 
bility, representing our ability to determine the system’s 
parameters through appropriate input/output measure¬ 
ments. 

(hi) Steering complex systems to desired states or tra¬ 
jectories. The ultimate goal of control is to drive a com¬ 
plex system from its current state/trajectory to some 
desired final state/trajectory. This problem has appli¬ 
cations from ecosystem management, to cell reprogram¬ 
ming. For example, we would like to design interven¬ 
tions that can move a cell from a disease (undesired) to 
a healthy (desired) state. We discuss different ways of 
achieving such control: (a) By applying small pertur¬ 
bations to a set of physically or experimentally feasible 
parameters; (b) Via compensatory perturbations of state 
variables that exploit the basin of attraction of the de¬ 
sired final state; (c) By mapping the control problem into 
a combinatorial optimization problem on the underlying 
network. 

(iv) Controlling collective behavior. Collective behav¬ 
ior, a much-studied topic in modern statistical physics, 
can result from the coordinated local activity of many 
interdependent components. Examples include the emer¬ 
gence of flocking in mobile agents or synchronization in 
coupled oscillators. Controlling such processes has nu¬ 
merous potential applications, from the design of flocking 


dent algorithm to find the minimum set of driver nodes robots ( 

Olfati-Saber 

2006 

, to the treatment of Parkin- 

required to control the network, nor an analytical frame- son’s disease ( 

Tass et al. 

1998 

). We review a broad 


work to estimate the fraction of driver nodes. Advances 
on this front became possible by mapping the control 
problem into well-studied network problems, like match¬ 
ing, and utilizing the notion of thermodynamic limit in 
statistical physics and the cavity method developed in 
spin glass theory, tools that were traditionally beyond 


the scope of control theory (Liu et al. 2011a). 


spectrum of methods to determine the conditions for the 
emergence of collective behavior and discuss pinning con¬ 
trol as an effective control strategy. 

Control problems are ubiquitous, with direct relevance 
to many natural, social and technological phenomena. 
Hence the advances reviewed here truly probe our fun¬ 
damental understanding of the complexity of the world 
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surrounding us, potentially inspiring advances in numer¬ 
ous disciplines. Consequently, our focus here is on con¬ 
ceptual advances and tools pertaining to control, that 
apply to a wide range of problems emerging in physical, 
technological, biological and social systems. It is this 
diversity of applications that makes control increasingly 
unavoidable in most disciplines. 


II. CONTROLLABILITY OF LINEAR SYSTEMS 


A system is controllable if we can drive it from any ini¬ 


tial state to any desired hnal state in hnite time (Kalman 


1963 1. Many mechanical problems can be formalized as 
controllability problems (Fig.[^. Consider, for example, 
the control of a rocket thrust upward. The rocket is con¬ 
trollable if we can find a continuous control input (thrust 
force) that can move the rocket from a given initial state 
(altitute and velocity) to a desired final state. Another 
example is the balancing of a stick on our hand. We know 
from our experience that this is possible, suggesting that 


the system must be controllable (Luenberger 19791. The 


scientific challenge is to decide for an arbitrary dynamical 
system if it is controllable or not, given a set of inputs. 

The current interest in the control of complex net¬ 
worked systems was induced by recent advances in the 


controllability of complex networks (Gao et al. 2014b 


Jia et aL\ |2013[ |Liu et aL\ 2011a 2012b Posfai et al. 


20131, offering mathematical tools to identify the driver 
nodes, a subset of nodes whose direct control with appro¬ 
priate signals can control the state of the full system. In 
general controllability is a prerequiste of control, hence 
understanding the topological factors of the underlying 
network that determine a system’s controllability offers 
numerous insights into the control principles of complex 
networked systems. As we discuss below, thanks to a con¬ 
vergence of tools from control theory, network science and 
statistical physics, our understanding of network control¬ 
lability has advanced considerably recently. 


A. Linear Time-Invariant Systems 

The starting point of most control theoretical ap¬ 
proaches is the linear time-invariant (LTI) control system 
(A,B) 


x(t) = Ax(t)-f B u(t). (11) 

Many mechanical systems can be naturally described 
by LTI dynamics, where the state vector captures the 
position and velocity of objects and the LTI dynamics is 
either directly derived from Newton’s Second Law or rep¬ 
resents some reasonable linearization of the underlying 
nonlinear problem, as illustrated by the stick balancing 
problem ([^. 


A = 


0 0 0 


021 0 0 
031 
041 
0 


0 0 
0 0 
0 0 
0 0 034 0 

0 055 

X — Ax -t- Bu 


0 0 


0 0 


; B = 



FIG. 3 (Color online) Graphical representation of a linear 
time-invariant system ( 111 . The state matrix A represents 


the weighted wiring diagram of the network that describes 
which components interact with each other and the direc¬ 
tion of the signal or information flow for each link; the input 
matrix B identifies the nodes (state variables) that are con¬ 
trolled by an outside controller. The network shown in the 
figure is controlled by an input vector u = (ui(t), U 2 (t))^ with 
two independent signals ui{t) and U 2 {t). The three actuator 
nodes {xi,X 2 and X 5 ) are the nodes directly controlled by 
u(t). These actuator nodes correspond to the three non-zero 
elements in B. The two driver nodes {xi and X2), represent¬ 
ing nodes that do not share input signals, correspond to the 
two columns of B. Note that node x^, is an actuator node, 
but not a driver node. 


A significant fraction of the control theory literature 
deals exclusively with linear systems. There are multiple 
reasons for this. First, linear systems offer an accurate 
model for some real problems, like consensus or agree¬ 
ment formation in multi-agent networks, where the state 


of each agent captures its opinion ( 

Liu et ai 

. 2008 

Mes- 

bahi and Egerstedtl |20101 |Rahmani et al. 2 

0 

0 

Tanner 

2004). Second, while many complex systems are char- 


acterized by nonlinear interactions between the compo¬ 
nents, the first step in any control challenge is to establish 


the controllability of the locally linearized system (Slo- 


tine and Li 1991). Furthermore, as we show below, for 


systems near their equilibrium points the linearized dy¬ 
namics can actually characterize the underlying nonlin¬ 
ear controllability problem. Third, the non-trivial net¬ 
work topology of real-world complex systems brings a 
new layer of complexity to controllability. Before we can 
explore the fully nonlinear dynamical setting, which is 
mathematically much harder, we must understand the 
impact of the topological characteristics on linear con¬ 
trollability, serving as a prerequisite of nonlinear control¬ 
lability. 

Consider the LTI dynamics a on a directed weighted 
network G(A) of N nodes (Fig. |^. The state variable 
Xi{t) can denote the amount of traffic that passes through 


a node i on a communication network (Pastor-Satorras 


and Vespignani 2004), or transcription factor concen¬ 


tration in a gene regulatory network (Lezon et al. 20061. 


The state matrix A := {aij)NxN represents the weighted 
wiring diagram of the underlying network, where is 
the strength or weight with which node j affects/influ¬ 
ences node i: a positive (or negative) means the link 
(j —>■ i) is excitatory (or inhibitory), and a^- = 0 if node 
j has no direct influence on node i. Consider M inde- 


































































7 


pendent control signals {mi, • • • , um} applied to the net¬ 
work. The input matrix B := {bim)NxM identifies the 
nodes that are directly controlled, where bim represents 
the strength of an external control signal Umif) injected 
into node i. 

The input signal u(t) = (ui{t), ■ ■ ■ ,UM{t))^ S can 
be imposed on all nodes or only a preselected subset of 
the nodes. In general the same signal Um{t) can drive 
multiple nodes. The nodes directly controlled by u(t) 
are called actuator nodes or simply actuators, like nodes 
xi, X 2 and x^ in Fig. The number of actuators is given 
by the number of non-zero elements in B. The actuators 
that do not share input signals, e.g. nodes Xi and X 2 
in Fig. are called driver nodes or simply drivers. The 
number of driver nodes equals the number of columns in 
B. 

Controllability, the ability to steer a system into an 
arbitrary final state in a finite time, implies that we can 
move the state variable of each node of a network to a 
predefined value, corresponding to the system’s desired 
position in the state space. Our ability to do so is greatly 
determined by the network topology. For example, if 
the network structure is such that a signal cannot get 
from our driver nodes to a particular node, that node, 
and hence the system as a whole, is uncontrollable. Our 
challenge is to decide when control is possible and when is 
not. The answer is given by controllability tests described 
next. 


B. Kalman's Criterion of Controllability 


Controllability tests allow us to check if an LTI system 
is controllable from a given set of inputs. The best known 


is Kalman’s rank condition (Kalman 1963), stating that 


the LTI system (A, B) is controllable if and only if the 
N X NM controllability matrix 


C = [B,AB,A^B,..,A^-^B] 
has full rank, i.e. 


( 12 ) 


rankC = N. (13) 


To understand the origin of ( |l^ , we consider the for¬ 
mal solution of ( 11 ) with x( 0 ) = 0 , i.e. 


x(t) = f exp[A(t — r)] B u(t) dr. 

Jo 


(14) 


If we expand exp[A(t — r)] in series, we will realize 
that x(t) is actually a linear combination of the columns 
in the matrices {B,AB,A2B, -}. Note that for any 
A' > N, we have rank [B, AB, A^B, • • • ,A^'-iB] = 
rankC. So if rankC < N, then even the infinite series of 
{B,AB,A2B, } will not contain a full basis to span 
the entire A-dimensional state space. In other words, we 
cannot fully explore the state space, regardless of u(t), 


indicating that given our inputs the system is stuck in a 
particular subspace, unable to reach an arbitrary point 
in the state space (Fig.|^. If, however, rankC = N, then 
we can find an appropriate input vector u{t) to steer the 
system from x(0) to an arbitrary x(t). Hence, the system 
is controllable. 

One can check that in the stick balancing problem ( [ 6 a| ) , 
the controllability matrix has full rank (rankC = N = 2), 
indicating that both systems are controllable. In the net¬ 
work control problem of Fig. & the controllability matrix 


C = 


6i 0 0 

0 02161 0 

0 03161 0 


(15) 


is always rank deficient, as long as the parameters 61 , 
021 E^nd O 31 are non-zero. Hence, the system is uncon¬ 
trollable. In contrast, for Fig. & we have 


C = 


61 0 0 0 0 0 

0 62 02161 0 0 0 

0 0 03161 0 0 0 


(16) 


which has full rank, as long as the parameters 61 , 62 , 021 
and O 31 are non-zero. Hence the system is controllable. 

The example of Fig. [^implies that the topology of the 
controlled network, which consists of both the network 
itself and the control signals applied to some nodes, im¬ 
poses some inherent limits on the controllability matrix: 
some configurations are controllable (Fig. [^), while oth¬ 
ers are not (Fig. |^). Thanks to the Kalman criterion, 
controllability can be easily tested when the dimension 
of the controllability matrix is small and its rank test can 
be done even without knowing the detailed values of its 
non-zero matrix elements. For large real networks the 
controllability test (131 is difficult to perform, however. 
Indeed, there is no scalable algorithm to numerically de¬ 
termine the rank of the controllability matrix C, which 
has dimension N x NM. Equally important, execut¬ 
ing an accurate rank test is ill-conditioned and is very 
sensitive to roundoff errors and uncertainties in the ma¬ 
trix elements. Indeed, if we plug the numerical values 
of bi and into ( 12 ), we may obtain extremely large 
or small matrix elements, such as which for large 

N are rather sensitive to numeric precision. Hence, for 
large complex systems we need to determine the system’s 
controllability without numerically calculating the rank 
of the controllability matrix. As we discuss in the next 
section, this can be achieved in the context of structural 
control theory. 


C. Structural Controllability 

For many complex networks the system parameters 
(e.g. the elements in A) are not precisely known. In¬ 
deed, we are often unable to measure the weights of the 













FIG. 4 (Color online) Controlling star networks, (a) Con¬ 
trolling the central node of a directed star does not assure 
controllability of the whole network, as shown in ( 15 l. (b) 


Indeed, the system is stuck in the plane 031*2(1) = 021*3(1), 
hence no signal Ui{t) can make the system leave this plane 
and explore the whole state space. The reason is simple: if 
we change Mi(t), *2(1) and *3(1) always evolve in a correlated 
fashion, indicating that we are unable to control the two nodes 
independently of each other. Note that while the system is not 
controllable in the whole state space, it remains controllable 
within the plane. It is natural that ensuring controllability 
within a restricted subspace will require fewer driver nodes 


than ensuring controllability within the whole state space (Liu 
et al. | 2011 b Muller and Schuppert[ [ 2011 1. (c) To ensure 


controllability, we must inject an additional signal U2 to ei¬ 
ther *2 or *3, in which case, according to (161, the network 
becomes controllable. After ( |Liu et ol j |20lib| ). 




FIG. 5 (Color online) Controllability, structural controllabil¬ 
ity, and strong structural controllability, (a) A directed path 
can be controlled by controlling the starting node only. The 
controllability is independent of the detailed (non-zero) val¬ 
ues of 61, a2i, and 032, so the system is strongly structurally 
controllable, (b) A directed star can never be controlled by 
controlling the central hub (node *1) only, (c) This network 
obtained by adding a self-edge to the star shown in b, can 
be controlled by controlling *1 only. The controllability is 
independent of the detailed (non-zero) values of 61, a2i, <131, 
and 033, so the system is strongly structurally controllable, 
(d) This network is controllable for almost all weights com¬ 
binations. It will be uncontrollable only in some pathological 
cases, for example when the weights satisfy the constraint 
032021 = 0230I1 exactly. Hence, the system is structurally 
controllable but does not display strong structural controlla¬ 
bility. 


links, knowing only whether there is a link or not. In 
other cases the links are time dependent, like the traf¬ 
fic on an internet cable or the flux of a chemical reac¬ 
tion. Hence, it is hard, if not conceptually impossible, to 
numerically verify Kalman’s rank condition using fixed 
weights. Structural control, introduced by C.-T. Lin in 
1970s, offers a framework to systematically avoid this 


limitation (Lin, 1974). 


1. The power of structural controllability 

An LTI system (A, B) is a structured system if the 
elements in A and B are either fixed zeros or indepen¬ 
dent free parameters. The corresponding matrices A and 
B are called structured matrices. The system (A, B) is 
structurally controllable if we can set the nonzero ele¬ 
ments in A and B such that the resulting system is con¬ 
trollable in the usual sense (i.e., rankC = N). 

The power of structural controllability comes from the 
fact that if a system is structurally controllable then it 
is controllable for almost all possible parameter realiza¬ 


tions (Davison, 1977 Dion et al.\ 2003 Glover and Silver- 
man[ 1976| Hosoe and Matsumoto[ |1979| |Lin[ |1974[ |Li 


nemann[ 1986| Mayeda[ 1981 Reinschke) 1988[ Shields 
and Pearson 19761. To see this, denote with S the set of 


all possible LTI systems that share the same zero-nonzero 
connectivity pattern as a structurally controllable sys¬ 
tem (A, B). It has been shown that almost all sys¬ 
tems that belong to the set S are controllable except for 


some pathological cases with Lebesgue measure zero (Lin 


1974 Shields and Pearsonf |1976 ). This is rooted in the 
fact that if a system (Aq, Bq) S 5 is uncontrollable, then 
for every e > 0 there exists a controllable system (A, B) 
with ||A — AqII < e and ||B — Bo|| < e where || • || de¬ 
notes matrix norm (Lee and Markus 1968 
In other words, an uncontrollable system in S becomes 
controllable if we slightly alter some of the link weights. 
For example, the system shown in Fig. [5ji is controllable 
for almost all parameter realizations, except when the 
edge weights satisfy the constraint 032021 = 023023 . But 
these pathological cases can be easily avoided by slightly 
changing one of the edge weights, hence this system is 
structurally controllable. 

Taken together, structural control tells us that we can 
decide a network’s controllability even if we do not know 
the precise weight of each edge. All we have to make 
sure is that we have an accurate map of the system’s 
wiring diagram, i.e., know which components are linked 
and which are not. As we demonstrate in the coming sec¬ 
tion, this framework considerably expands the practical 
applicability of control tools to real systems. 


2. Graphical interpretation 

Structural control theory allows us to check if a con¬ 
trolled network is structurally controllable by simply in¬ 
specting its topology, avoiding expensive matrix oper¬ 
ations. This is possible thanks to the graphical inter- 
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pretation of Lin’s Structural Controllability Theorem, 
discussed next. 

Consider an LTI system (A, B) represented by a di¬ 
graph G'(A,B) = (y,E) (Fig. [^. The vertex set 
V = Va ^ Vb includes both the state vertices Va = 


{xi,--- jXat} = {ui,--- ,vn}, corresponding to the N 
nodes of the network, and the input vertices Vb = 
{ui,--- ,um} = {uAT-i-i,-- - ,vn+m}, corresponding to 
the M input signals that are called the origins or roots 
of the digraph G(A,B). The edge set E = Ea U Eb 
includes both the edges among state vertices Ea = 
{{xj,Xi)\aij 7 ^ 0}, corresponding to the links of network 
A, and the edges connecting input vertices to state ver¬ 
tices Eb = {ium,Xi)\bim 7 ^ 0}. These definitions allow 
us to formulate a useful statement: The system (A, B) 
is not structurally controllable if and only if it has inac¬ 
cessible nodes or dilations (Lin 19741. 

Let us consider these two cases separately. A state ver¬ 
tex Xi is inaccessible if there are no directed paths reach¬ 
ing Xi from the input vertices (Fig. Consequently, an 
inaccessible node can not be influenced by input signals 
applied to the driver nodes, making the whole network 
uncontrollable. 

The digraph G(A, B) contains a dilation if there is a 
subset of nodes S C Va such that the neighborhood set 
of S, denoted as T{S), has fewer nodes than S itself (see 
Fig. 1^). Here, T{S) is the set of vertices Vj for which 
there is a directed edge from Vj to some other vertex in S. 
Note that the input vertices are not allowed to belong to 
S but may belong to T{S). Roughly speaking, dilations 
are subgraphs in which a small subset of nodes attempts 
to rule a larger subset of nodes. In other words, there are 
more “subordinates” than “superiors”. A controlled net¬ 
work containing dilations is uncontrollable. For example, 
in a directed star configuration, where we wish to control 
via a central node all the leaves, any two leaf-nodes form 
a dilation with the central hub. If we control the central 
hub only, the system remains uncontrollable because we 
cannot independently control the difference between the 
two leaf nodes’ states (Fig.|^. In other words, we cannot 
independently control two subordinates if they share the 
same superior. 

Taken together, Lin’s structural controllability theo¬ 
rem states that an LTI system (A, B) is structurally con¬ 
trollable if and only if the digraph G(A, B) does not con¬ 
tain inaccessible nodes or dilations. These two conditions 
can be accurately checked by inspecting the topology 
of the digraph G(A, B) without dealing with floating¬ 
point operations. Hence, this bypasses the numerical is¬ 
sues involved in evaluating Kalman’s controllability rank 


^ The structural controllability theorem also has a pure algebraic 
meaning ( [Shields and Pearson[ [1976^ , which plays an impor¬ 
tant role in the characterization of strong structural controlla¬ 
bility (jMayeda and Yamadaj [l979|. 



FIG. 6 (Color online) Inaccessibility, dilations and cacti, (a) 
The red nodes {xi,X2) are inaccessible from the input node 
ui (in blue), as variations in ui do not influence the state of 
Xi or X2. (b) The red nodes S = {x3,Xii} cause a dilation. 
Indeed, their neighborhood set T{S) = {*5} contains only 
one node, hence the size of T{S) is smaller than S, implying 
that a single node in T( 5 ) aims to control two nodes in S. As 
we showed in Eq. ( 151 and Fig. this is not possible, (c) A 
cactus contains neither inaccessible nodes nor dilations. Note 
that in the cactus structure T{S) = {x2,X5}, hence there is 
no dilation. There is only one stem (shown in green) in one 
cactus. There could be multiple buds (shown in purple) in the 
same cactus. A cactus is a minimal structure for structural 
controllability. 


test, and also our lack of detailed knowledge on the edge 
weights in G(A,B). 

An alternative graph theoretical formulation of Lin’s 
structural controllability theorem is often useful in prac¬ 
tice. A general graph is covered or spanned by a subgraph 
if the subgraph and the graph have the same vertex set. 
Typically the spanning subgraph has only a subset of 
links of the original graph. For a digraph, a sequence 
of oriented edges {(ui —)■ ^ 2 ), - • • , (vk-i —f Vk)}, where 
the vertices {ui,U 2 , • • • ,Vk} are distinct, is called an ele¬ 
mentary path. When Vk coincides with vi, the sequence 
of edges is called an elementary cycle. For the digraph 
G(A, B), we define the following subgraphs (Fig.|^): (i) 
a stem is an elementary path originating from an input 
vertex; (ii) a bud is an elementary cycle G with an addi¬ 
tional edge e that ends, but does not begin, in a vertex 
of the cycle; (iii) a cactus is defined recursively: A stem 
is a cactus. Let G, O, and e be, respectively, a cactus, 
an elementary cycle that is disjoint with G, and an arc 
that connects G to O in G(A, B). Then, G U {e} U O is 
also a cactus. G(A, B) is spanned by cacti if there exists 
a set of disjoint cacti that cover all state vertices. 


Note that a cactus is a minimal structure that con¬ 
tains neither inaccessible nodes nor dilations. That is, 
for a given cactus, the removal of any edge will result in 
either inaccessibility or dilation, hence the controllabil¬ 
ity of the cactus is lost (Fig. |^. We can now formulate 
Lin’s structural controllability theorem as follows: An LTI 
system (A, B) is structurally controllable if and only if 
G(A,B) is spanned by cacti ( |Lin[ |I974[ ) . Later we show 
that this formulation helps us design an efficient algo¬ 
rithm to identify a minimum set of inputs that guarantee 
structural controllability. 
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3. Strong structural controllability 


The fundamental assumption of structural control is 
that the entries of the matrices A and B are either ze¬ 
ros or independent free parameters. Therefore struc¬ 
tural control does not require knowledge of the exact 
values of parameters, any by avoiding floating-point op¬ 
erations, it is not subject to numerical errors. However, 
some systems have interdependent parameters, making 
it uncontrollable despite the fact that it is structurally 
controllable. For example. Fig. displays an LTI sys¬ 
tem that is structurally controllable, but becomes un¬ 
controllable when the parameters satisfy the constraint 
0320.21 = a 23 a§i. This leads to the notion of strong 
structural controllability (SSC) : A system is strongly 
structurally controllable if it remains controllable for any 
value (other than zero) of the indeterminate parame¬ 
ters ( Mayeda and Yamada[ 1979). In other words, there 
is no combination of non-zero link weights that violates 


Kalman’s criterion (13). For example, the LTI systems 


shown in Fig. and c are strongly structurally control¬ 
lable. 


and Yamada 


schke et al. 


Both graph-theoretic (Jarczyk et al. 2011 Mayeda 


19791 and algebraic conditions (Rein 


1992) for SSC have been studied. Un¬ 


fortunately, those conditions do not lead to efficient 
algorithms. Recently, necessary and sufficient graph- 
theoretical conditions involving constrained matchings 


were derived (Chapman and Mesbahi 2013). Denote a 


matching of size t in the bipartite representation H{A) 
of the digraph G(A) as t-matching. A t-matching is con¬ 
strained if it is the only t-matching in H(A). A match¬ 
ing is called I4-less if it contains no edges correspond¬ 
ing to self-loops. Let S be an input set with cardinality 
M < N. The corresponding structured pair (A,B) is 
strongly structurally controllable if and only if H{A) has 
a constrained {N — M)-matching with S unmatched and 
i7(Ax) has a constrained 14-less (A —M)-matching with 
S unmatched. Here H{A^) is formed by adding self¬ 
loops to all nodes if they don’t have one. The constrained 
matching conditions can be applied to check if an input 
set is strongly structural controllable in 0{N‘^). Though 
finding a minimum cardinality input set is proven to be 
NP-complete, a greedy 0{N^) algorithm has been de¬ 
veloped to provide a strongly structural controllable in¬ 


put set, which is not-necessarily minimal (Chapman and 


Mesbahi 2013). 


D. Minimum Input Problem 

If we wish to control a networked system, we first need 
to identify the set of driver nodes that, if driven by differ¬ 
ent signals, can offer full control over the network. Any 
system is fully controllable if we control each node indi¬ 
vidually. Yet, such full control is costly and typically im¬ 


practical. Hence, we are particularly interested in identi¬ 
fying a minimum driver node set (MDNS), whose control 
is sufficient to make the whole system controllable. In 
other words, we want to control a system with minimal 
inputs. 


1. Solution based on structural control theory 


Kalman’s rank condition does not offer us the MDNS— 
it only tells us if we can control a system through a 
given set of potential driver nodes that we must guess 
or select. Furthermore, to numerically check Kalman’s 
rank condition, we have to know all the entries in A 
and B, which are often unknown for complex networks. 
Even if we know all the weights (parameters) exactly, 
a brute-force search for the MDNS would require us 
to compute the rank of almost 2^ distinct controllabil¬ 
ity matrices, a combinatorially prohibitive task for any 
network of reasonable size. Yet, as we show next, we 
can identify the MDNS by mapping the control prob¬ 
lem into a purely graph theoretical problem called maxi- 


mum matching (Commault et al. 

o 

o 

Liu et al. 

2011 a 

Murota 

2009 Yamada and Fould 

s 199 

0 ). 


Matching is a widely studied problem in graph theory. 


with many practical applications (Lovasz and Plummer 


2009). On undirected graphs, where it was originally de¬ 


fined, a matching represents a set of edges without com¬ 
mon vertices (red edges in Fig. [^). Maximum match¬ 
ing is a matching of the largest size. For most graphs 
we can find multiple maximum matchings (Fig. -h3). 
The end vertices of a matching edge are called matched, 
the remaining vertices are unmatched. If all vertices are 
matched, then the matching is perfect (Fig. B)- 

Many real world problems can be formalized as a max¬ 
imum matching problem on bipartite graphs (Fig. 
Consider, for example, M job applicants applying for N 
openings. Each applicant is interested in a subset of the 
openings. Each opening can only accept one applicant 
and an applicant can only accept one job offer. Find¬ 
ing an assignment of openings to applicants such that 
as many applicants as possible get a job is a classical 
maximum matching problem. 

In structural control theory, the role of matching is 
well studied and matching was originally defined in the 


bipartite representation of a digraph ( 

Commault et al. 

2002 

Murota 

2009 Yamada and Foulds 

19901. The ex- 

tended definition of matching on a digraph (Liu et al. 


2011 a) connects more naturally to the cactus structure 
(Fig. which is a fundamental notion in structural con¬ 
trol theory. In a directed graph (digraph), a matching 
is defined to be a set of directed edges that do not share 


common start or end vertices (Liu et al.. 2011a). Hence, 


a vertex can be the starting or the end point of a red 
link, but we cannot have two red links pointing to the 
same vertex. A vertex is matched if it is the end ver- 
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FIG. 7 (Color online) Matching. The figures show the maximum matchings of (a,b) undirected graphs, (c) a bipartite graph 
and (d,e,f) digraphs. For undirected or bipartite graphs, a matching represents a set of edges without common vertices. For 
digraphs, a matching is a set of directed edges that do not share the common start or end vertices. Maximum matching is a 
matching with the largest number of edges. On panels (g-i) edges in the matching are colored in red. Matched (or unmatched) 
nodes are shown in green (or white), respectively. 


tex of a matching edge. Otherwise, it is unmatched. For 
example, in a directed path, all but the starting vertex 
are matched (Fig. ,j). A matching of maximum size 
is called a maximum matching. A maximum matching 
is called perfect if all vertices are matched, like in a di¬ 
rected elementary cycle (Fig. [^,1). We can prove that 
a matching of a digraph can be decomposed into a set 
of directed paths and/or directed cycles (Fig. §r). Note 
that directed paths and cycles are also the basic elements 
of the cactus structure (Fig. |^). Hence, matching in di¬ 
graphs connects naturally to the cactus structure. 

The usefulness of matching in network control comes 
from a theorem that provides the minimum number of 


driver nodes in a network (Liu et al. 2011a). 


Minimum input theorem: To fully control a directed 
network G(A), the minimum number of inputs, or equiv¬ 
alently the minimum number of driver nodes, is 


Ad = max{iV- |M*|,1}, 


(17) 


where \M* \ is the size of the maximum matching in G(A). 
In other words, the driver nodes correspond to the un¬ 
matched nodes. If all nodes are matched {\M*\ = N), 
we need at least one input to control the network, hence 
Ad = 1. We can choose any node as our driver node in 
this case. 

The minimum input theorem maps an inherently dy¬ 
namical problem, i.e. our ability to control a network 
from a given subset of nodes, into a purely graph the¬ 
oretical problem of finding the maximum matching of a 
directed network. Most important, it bypasses the need 
to search all node combinations for a minimum driver 
node set, as the driver nodes are provided by the solu¬ 
tion of the underlying matching problem. 

Maximum matching: algorithmic solution. The map¬ 
ping of the MDNS problem to a matching problem via 


(17) seems to map a problem of high computational corn- 


maximum matching for a digraph. The real value of this 
mapping, however, comes from the fact that the max¬ 
imum matching problem in a digraph in not NP-hard, 
but can be solved in polynomial time. Indeed, the max¬ 
imum matching for a digraph can be identified by map¬ 
ping the digraph to its bipartite representation, as illus¬ 
trated in Fig. Consider a digraph G(A), whose bi¬ 
partite representation is H{A) = {V^ U V^,T). Here, 
Vj' = {cc/", • • • , x/)} and VJ[ = {xf, • • • , are the 
set of vertices corresponding to the A columns and rows 
of the state matrix A, respectively. The edge set of this 
bipartite graph is F = {{x^,x~) \ aij ^ 0}. In other 
words, we split each node Xi of the original digraph into 
two “nodes” xf and x~. We then place an edge (x^, x~) 
in the bipartite graph if there is a directed edge (xj Xi) 
in the original digraph. Note that since we allow self¬ 
loops (xi —)■ Xi) in the original digraph, there can be 
edges of this type (x/',x“) in the bipartite graph. A 
maximum matching of a bipartite graph can be found ef¬ 
ficiently using the Hopcroft-Karp algorithm , which runs 
in 0{'/VE) time (Hopcroft and Karp, 1973). After run¬ 
ning the algorithm, we can map the maximum matching 
in the bipartite representation, e.g. {x^ ,xf), (x^, Xg ) in 
Fig.[9]o, back to the maximum matching in the original 
diagraph, e.g. (xi, X 2 ), (xg, Xg) in Fig. [^, obtaining the 
desired maximum matching and hence the corresponding 
MDNS. 

Taken together, the maximum matching algorithm al¬ 
lows the efficient identification of the MDNS using the 


plexity — an exhaustive search for the MDNS — into 
another just as complicated problem, that of finding the 


following steps (Fig. 10): (i) Start from the directed net¬ 
work we wish to control and generate its bipartite repre¬ 
sentation (Fig.j^. Next identify a maximum matching on 
the underlying bipartite graph using the Hopcroft-Karp 
algorithm, (ii) To each unmatched node add a unique 
control signal, as unmatched nodes represent the driver 
nodes, (iii) As there could be multiple maximum match¬ 
ings for a general digraph, multiple MDNSs exist, with 
the same size Ad. 
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Maximum Matching U-rooted factorial connection Cacti Structure Controlled Network 


FIG. 8 (Color online) Graph-theoretic proof of the minimum input theorem, (a) A directed network, (b) The maximum 
matching represents the largest set of edges without common heads or tails. All maximum matchings can be decomposed into a 
set of vertex-disjoint directed paths and directed cycles, shown in red. If a node is the head of a matching edge, then this node 
is matched (shown in green). Otherwise, it is unmatched (shown in white). The unmatched nodes must be directly controlled 
to control the whole network, hence they are the driver nodes, (c) By injecting signals into driver nodes, we get a set of directed 
paths whose starting points are the input nodes. The resulting paths are called “stems” and the resulting digraph is called 
U-rooted factorial connection, (d) By “grafting” the directed cycles to those “stems”, we get “buds”. The resulting digraph 
is called cactus or cacti. A cactus is a minimal structure for structural controllability, as removing any of its edges will cause 
either inaccessible nodes or dilations, (e) According to the structural controllability theorem, since there is a cacti structure 
(highlighted in yellow) underlying the controlled network, the system is structurally controllable. Note that (a-d) also suggests 
an efficient method to identify the minimal cacti, i.e. the cacti structure with the minimum number of roots. This minimal 
cacti serve as the control skeleton that maintains the structural controllability of the system. 



b 



FIG. 9 (Golor online) Maximum matching calculation. The 
maximum matching of the digraph (a) can be computed from 
its bipartite representation (b), which is obtained by splitting 
each node Xi into two “nodes” {xf and x~) and placing an 
edge {x'j ,x~) in the bipartite graph if there is a directed edge 
{xj —^ Xi) in the original digraph. The maximum matching of 
any bipartite graph can be identified in polynomial time using 
the Hopcroft-Karp algorithm. Mapped back to the digraph, 
we obtain the maximum matching of the original digraph and 
the driver nodes of the corresponding control problem. 



FIG. 10 (Golor online) Identifying the driver nodes. For a 
general directed network, like the one shown in the left panel, 
there could be multiple maximum matchings, shown in red 
on the right panels. Hence, we can identify multiple MDNSs 
(white nodes). To each driver node we must add a unique 
control signal necessary to ensure structural controllability. 


Recently, several algorithmic approaches have been de¬ 
veloped to optimize the network controllability (in the 


sense of decreasing Nd) via minimal structural pertur¬ 
bations, like adding a minimum number of edges at ju¬ 


diciously chosen locations in the network (Wang et al. 


2012), rewiring redundant edges (Hou et al.\ |2013|, and 


assigning the direction of edges (Hou et aZ.l 12012 


et al. , 2014). 


Xiao 


Maximum matching: analytical solution based on the 
cavity method. While the maximum matching allows 
us to efficiently identify the MDNS, the algorithmic ap¬ 
proach provides no physical insights about the impact 
of the network topology on N^. For example, what 
network characteristics influence N^, and how does Nd 
depend on them? Which networks are easier to con¬ 
trol and which are harder? To answer these questions 
we can turn to the cavity method, a versatile tool of 


statistical physics (Mezard and Parisi, 2001 Zdeborova 
and Mezard 2006 Zhou and Ou-Yang 2003[ ). We il- 
lustrate this approach by analytically calculating np), 
representing the fraction of driver nodes ud (= Nd/N) 
averaged over all network realizations compatible with 
the network’s degree distribution P{kin, fcout) ( |Liu et al] 
2011a). We start by describing a matching M in a 


digraph G = {V{G), E{G)} by the binary variables 
Sa = ^ {0; 1} assigned to each directed edge 

a = {i ^ j) G E[G) with = 1 if a belongs to the 
matching M and Sa = 0 otherwise. According to the def¬ 
inition of matching in a digraph, matching edges do not 
share starting or end nodes, formally resulting in two con¬ 
straints for each vertex i G V(G): (i) ^(^->- 3 ) — 

(ii) SfcGO-i ^ with d~i and d~^i indicating the 

sets of nodes that point to i or are pointed by i, respec¬ 
tively. 


The quantity fi({s}) = 1 — J2k&d-i tells us the 

state of each vertex: vertex i is matched if fids}) = 0 
and unmatched if fi({s}) = 1. Consequently, the cost 
(or energy) function gives for each matching M = {s} 
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the number of unmatched vertices 

£g({4)= E ms})=N-\M\. (18) 

ieV{G) 

We define the Boltzmann probability in the space of 
matchings as 


Vg{{s}) = 


ZgW) ' 


(19) 


where j3 is the inverse temperature and Zq{P) is the par¬ 
tition function 


zgW) = J2‘ 


-/3fG({s}) 


( 20 ) 


{«} 


In the limit /3 —>■ oo (i.e. the zero temperature limit), 
the internal energy SgW) and the entropy SciP) provide 
the ground state properties, i.e. the properties of the 
maximum matchings. In particular, Scioo) represents 
the number of unmatched vertices (with respect to any 
maximum matching), and the entropy Sg{oo) yields the 
logarithm of the number of maximum matchings. 

In the zero temperature limit, the average fraction of 
driver nodes is given by 

no = [G(ty2) + 0(1 - wi) - l] 

+ [G(n;2) + G(l-n;i)-l] 

+ 1 [wi(l - 1 (^ 2 ) + wi(l - W 2 )] |, ( 21 ) 

where Wi,W 2 ,W 3 ,Wi,W 2 ,W 3 satisfy the set of self- 
consistent equations 


Wi = H{w2) 

W2 = I — H{1 — Wi) 
W 3 = 1 — W 2 — Wi 

Wi = H{w2) 

W 2 = ^ — H[1 — rci) 
W 3 = 1 — W 2 — Wl 


( 22 ) 


and 


'G(x) = 

G(a:) = EZ=o Pihn)x'^'’^ 

= Sfcout=o Qi^ont + l)a;^°"* 

.H{x)= EZ=oQi^in + 


(23) 


are the generating functions, and Q{kont) = ^ 

Q(fci„) = are the out- and in- degree distribu¬ 

tions of the node i when one selects uniformly at random 
a directed edge {i ^ j) from the digraph. 

While the cavity method does not offer a closed-form 


impact of key network characteristics, like the average 
degree (fc) or the degree exponent 7 , on ud in the ther¬ 
modynamic limit {N —>■ 00 ). For example, for directed 
Erdos-Renyi random networks ( [Bollobas 2001 Erdos 
and Renyi, 1960), both P{kin) and P(A:out) follow a Pois¬ 
son distribution, i.e. e~''^'>^‘^{{k)/2Y/k\. In the large {k) 
limit we have 


riD 




(24) 


For directed scale-free networks, we assume that 
P{ki-a) and P{koE) have the same functional form with 
power-law exponent 7 and exponential cutoff P{ki^) = 
GfcWe”^/'^, P(fcout) = C . Here the normal¬ 
ization constant is G = [L^(e“^/'')] , where Li„(a;) is 

the nth polylogarithm of x. Due to the exponential cut¬ 
off the distribution is normalizable for any 7 . One 

can show that as 7 —>■ 2, we have ud —>■ 1. This means 
one has to control almost all the nodes to achieve full 
control over the network. Therefore 7 = 2 is the critical 
value for the controllability of scale-free networks, as only 
for 7 > 2 can we obtain full controllability by controlling 
only a subset of the nodes. Note that for 7 —)■ 2 super¬ 
hubs emerge that connect to almost all nodes in the net¬ 
work (Albert and Barabasi 2002 Barabasi 2015). We 


know that for a star-like digraph with one central hub 
and — 1 leaves, one has to control A^d = A^ — 1 nodes 
(the central hub and any N — 2 leaves). In the large N 
limit, A^ — 1 « A^, which explains intuitively why we have 
to control almost all nodes when 7 —?> 2 . 

For scale-free networks with degree exponent 7 in = 


7 out = 7 generated from the static model (Goh et al. 


2001), the parameters (fc) and 7 are independent. In the 


thermodynamic limit the degree distribution is P{k) = 
^ where r(s) is the gamma 

function and r(s,x) the upper incomplete gamma func¬ 
tion. In the large k limit, P{k) ~ where 

7 = 1-1-^. The asymptotic behavior of nuUk),^) for 
large (fc) is 


no 




(25) 


If 7in 7^ 7out; the smaller of the two exponents, i.e. 
min[ 7 in, 7out] determines the asymptotic behavior of nu- 
Equation (25) indicates that as 7 —)■ 2, nu —7 1, which 


is consistent with the result that 7 c = 2 for a purely SF 
network. 

The systematic dependence of ud on (fc) and 7 prompts 
us to ask: How do other network characteristics, like de¬ 
gree correlations, clustering, modularity, or the fraction 
of low degree nodes, influence nu (Menichetti et al. 


2014 Posfai et al. 2013). A combination of analytical 


solution, Eq. (21) allows us to systematically study the 


and numerical results indicate that the clustering coeffi¬ 
cient and modularity have no discernible effect on ud . At 
the same time the symmetries of the underlying match¬ 
ing problem generate linear, quadratic or no dependence 
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FIG. 11 (Color online) Analytical results on the fraction of 
driver nodes (wd = Nu/N) for canonical model networks, (a) 
For directed Erdds-Renyi random networks, no decays expo¬ 
nentially for large (k). (b) For directed scale-free networks no 
approaches one as the degree exponent 7 approaches two, in¬ 
dicating that in such networks all nodes need to be controlled. 


on degree correlation coefficients, depending on the na¬ 
ture of the underlying degree correlations (]P6sfai et al. 


20131. 


For uncorrelated directed networks, the density of 
nodes with fcin, fcout = 1 or 2 determine the size of maxi¬ 
mum matchings (Menichetti et aZ. 2014). This suggests 
that random networks whose minimum and Zcout are 
greater than two typically have perfect matchings and 
hence can be fully controlled via a single control input 
(i.e. Nu = 1), regardless of the other properties of the 
degree distribution. 


2. Solution based on PBH controllability test 


In structural control theory we assume that the system 
parameters, like the link weights in G(A, B), are either 
fixed zeroes or independent free parameters. This frame¬ 
work is ideal for many systems for which we only know 
the underlying wiring diagram (i.e. zero/nonzero val¬ 
ues, indicating the absence/presence of physical connec¬ 
tions) but not the link characteristics, like their weights. 
Yet, the independent free parameter assumption is very 
strong, and it is violated in some systems, like in undi¬ 
rected networks, where the state matrix A is symmet¬ 
ric, or unweighted networks, where all link weights are 
the same. In such cases structural control theory could 
yield misleading results on the minimum number of driver 
nodes N^. Hence, it is important to move beyond struc¬ 
tural control as we explore the controllability and other 
control related issues. 

For LTI systems with exactly known system parame¬ 
ters the minimum input problem can be efficiently solved 
using the Popov-Belevitch-Hautus (PBH) controllability 
test. The PBH controllability test states that the system 
(A,B) is controllable if and only if (Hautus 1969|) 


rank [si - A, B] = A, Vs G C. 


(26) 


Since the first N x N block of the N x {N + M) matrix 
[si - A, B] has full rank whenever s is not an eigenvalue 


of A, we only need to check each eigenvalue of A, i.e. 
s G A (A), when running the PBH test. 

Note that the PBH test (261 and Kalman’s rank condi¬ 
tion ( |T^ are equivalent. Yet, the advantage of the PBH 
test comes from the fact that it connects the controllabil¬ 
ity of (A, B) to the eigenvalues and eigenvectors of the 
state matrix A. This can be used to solve the minimum 
input problem exactly. Indeed, the PBH controllability 
test suggests that (A, B) is controllable if and only if 
there is no left eigenvector of A orthogonal to all the 
columns of B. In other words, the columns of B must 
have a component in each eigendirection of A. Recall 
that for an eigenvalue Aq G A(A), its algebraic multiplic¬ 
ity is the multiplicity of Aq as a root of the characteristic 
polynomial p(A) = det(A—AI). Its geometric multiplicity 
is the maximal number of linearly independent eigenvec¬ 
tors corresponding to it. Hence, the number of control 
inputs must be greater than or equal to the largest ge- 


ometric multiplicity of the eigenvalues of A ( 

Antsaklis 

and Michel 

1997 

Sontag 

1998 

Yuan et al. 

2013). In 


other words, the minimum number of control inputs (or 
equivalently the minimum number of driver nodes) is de¬ 
termined by the maximum geometric multiplicity of the 
eigenvalues of A, i.e. 


Ad = max{^(Ai)}, (27) 

i 

where /r(Ai) = dimVAi = A — rank(AilAr — A) is the 
geometric multiplicity of A’s eigenvalue A^, representing 
the dimension of its eigenspace. Note that the algebraic 
multiplicity of eigenvalue Ai, denoted by i5(Ai), is its mul¬ 
tiplicity as a root of the characteristic polynomial. In 
general, 5{Xi) > /i(Ai). But for symmetric A, which is 
the case for undirected networks, we have 5{\i) = /J.(Ai). 


TABLE I Eigenvalues and minimum number of driver nodes 
of some special graphs of A nodes. For unweighted and undi¬ 
rected star and fully connected networks, the table shows the 
algebraic multiplicity of eigenvalues in the parenthesis. After 
(Yuan et al. 20131. 


Network 

Eigenvalue 


Ad 

Chain 

2cos^, 1,. 

•• ,A 

1 

Ring 

2cos5^,g = l 

,••• ,A 

2 

Star 

0(A- 2),±^A^ 

1(1) 

A-2 

Complete graph A — 1(1), —1(A — 

1) 

A- 1 


Based on (27), we can develop an efficient algorithm 
to identify the minimum set of driver nodes for arbitrary 
LTI systems (Fig. 12), allowing us to explore the impact 
of the network topology and link-weight distributions on 
Ad (Yuan et aZ. , |2013 ). For undirected and unweighted 
ER networks of connectivity probability p, the results in¬ 
dicate that for small p, nu decreases with p, while for 
sufficiently large p, ud increases to (A — 1)/A, which 
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elementary column Column Canonical 
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FIG. 12 (Color online) Identifying a minimum set of driver 
nodes of small networks. For each network, we show the ma¬ 
trix A — A^'^I, its column canonical form, all eigenvalues A of 
A, and the eigenvalue A^^ with the largest geometric multi¬ 
plicity. We highlight the rows that are linearly dependent on 
others in the column canonical form in red. The correspond¬ 
ing nodes are the driver nodes (shown in red) of the corre¬ 
sponding networks. For undirected networks in (a) and (c), 
/r(A^'^) is equal to the maximum algebraic multiplicity, that 
is, the multiplicity of A*^. The conhguration of driver nodes 
is not unique as it relies on the elementary column trans¬ 
formation, but the minimum number of drivers is uniquely 
determined by the maximum geometric multiplicity /i(A*^) of 
matrix A. After | |Yuan et al. 1 [M^ . 


is exact for p = I (complete graph, see Table. |T]). This 
approach has been recently extended to multiplex net¬ 


works (Yuan et al.\ 2014|. 


E. Minimal Controllability Problems 


Any networked system with LTI dynamics is fully con¬ 
trollable if we control each node individually with an in¬ 
dependent signal, i.e. M = N. But this is costly and 
typically impractical for large complex systems. Hence, 
we are particularly interested in fully controlling a net¬ 
work with minimum number of nodes. Depending on the 
objective function and the way we “inject” input signals, 
we can formalize different types of minimal controllahility 
problems (MCPs) ( [Olshevsky 2014). 

(MCPO): One scenario is that we try to minimize the 
number of independent control signals, corresponding to 
the number of columns in the input matrix B, or equiv¬ 


alently, the number of driver nodes (Liu et al. 201 la I 


whose control is sufficient to fully control the system’s 


dynamics (Fig. 13 1 ). This is nothing but the minimum 


inputs problem discussed in the previous subsection. 

(MCPl): We assume dedicated inputs, i.e. each con¬ 
trol input Ui can only directly control one node (state 
variable). In the matrix form, this amounts to hnding 
a diagonal matrix B G that has as few nonzero 


3 MCPO: Minimize the number of driver nodes. 



X o' 
0 X 
X 0 
X 0 
0 0 


[3 MCP1: Minimizer the number of actuator nodes. 

Ul U4 

U3 f U2 


,-#Xi ' ’ -^4V., 1 ’ 

/ 'h H 


C MCP2: Minimizer the number of actuator nodes. 



FIG. 13 (Color online) Different minimal controllability prob¬ 
lems (MCPs). For each MCP, we show the corresponding 
graph representation G(A, B), and the input matrix B (where 
X’s stand for non-zero elements. MCPO: We aim to minimize 
the number of driver nodes, or equivalently, the number of in¬ 
dependent input signals. One signal can drive multiple nodes. 
MCPl: We aim to minimizer the number of actuator nodes 
which receive input signals. One signal can only drive one 
actuator node. MCP2: We aim to minimizer the number of 
actuator nodes with only one signal. This unique signal can 
drive multiple actuator nodes. In all cases, we assume there 
are four actuator nodes (xi,X 2 ,X 3 and 3 : 4 ). We color the 
driver nodes in pink. 


entries as possible so that the LTI system x = Ax + Bu 
is controllable (Fig. [T^). 

(MCP2): We set Ui(t) = u{t) and aim to find a vector 
b that has as few nonzero entries as possible such that 
the system x = Ax -f bu is controllable (Fig. 13:). 


Note that in solving MCPO, one signal can be applied 
to multiple nodes. The number of actuator nodes (corre¬ 
sponding to those non-zero entries in B) is not necessarily 
minimized. In MCPl u(t) is a vector of control inputs, 
i.e. we have multiple input signals, while in MCP2, u{t) 
is a scalar, i.e. there is only one input signal. In both 
cases, we try to minimize the number of actuator nodes 
that are directly controlled by input signals. 

Though MCPO for a general LTI system is easy to 
solve, MCPl and MCP2 are NP-hard ( [Olshevsky 2014). 
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Yet, if we need to guarantee only structural controlla¬ 
bility, MCPl can be easily solved (Pequito et al. 

20161. For a directed network G with LTI dynamics 


the minimum number of dedicated inputs (or actuators), 
Nda, required to assure structural controllability, is 


.^da — .^D + /3 — a. 


(28) 


where Yd is the minimum number of driver nodes; /3 
is the number of root strongly connected components 
(rSCCs), which have no incoming links from other SCCs; 
and a is the maximum assignability index of the bipar¬ 
tite representation B{G) of the directed network G. An 
rSCC is said to be a top assignable SCC if it contains at 
least one driver node with respect to a particular max¬ 
imum matching M*. The maximum assignability index 
of B{G) is the maximum number of top assignable SCCs 
that a maximum matching M* may lead to. The mini¬ 
mum set of actuators can be found with polynomial time 


complexity (Pequito et al.. 2013 20161. 


Consider, for example, the network shown in Fig. 
which has two possible maximum matchings Mi = 

{(xi —>■ X4),(X4 —>■ X3),(x5 —>■ X5)}, M2 = {(xi —>■ 

X2), (x 4 — >■ X3), (x 5 — >■ X5)}. Both have size 3, hence the 
number of driver nodes is Yd = max{Y — |M*|, 1} = 2, 
according to ( [l7| . Note that the two maximum match¬ 
ings will yield two minimum sets of driver nodes, i.e. 
{xi,X 2 } and {xi,X 4 }. The former is shown in Fig. 
There are two rSCCs, {xi} and {xs}, each containing 
a single node, hence /3 = 2. The rSCC {xi} is a top 
assignable SCC, because it contains one driver node with 
respect to either Mi or M 2 . The rSCC {xs} is not a 
top assignable SCC, because it contains no driver nodes. 
Hence the maximum assignability index of this system 
is q; = 1. Finally, the minimum number of actuators is 
Ya = Yd -I- /3 — a = 3 and there are two minimum sets 
of actuators, i.e. {xi,X 2 ,X 5 } and {xi,X 4 ,X 5 }. 


F. Role of Individual Nodes and Links 


As we have seen in Sec lII.P.Tj a system with Yd driver 
nodes can be controlled by multiple driver node con¬ 
figurations, each corresponding to a different maximum 


matching (Fig. 10). Some links may appear more often 


in the maximum matchings than other links. This raises 
a fundamental question: What is the role of the individ¬ 
ual node (or link) in control? Are some nodes (or links) 
more important for control than others? To answer these 
questions, in this section we discuss the classification of 
nodes and links based on their role and importance in 
the control of a given network. 


1. Link classification 

In both natural and technological systems we need to 
quantify how robust is our ability to control a network 


under unavoidable link failure. To adress this question, 
we can use structural controllability to classify each link 
into one of the following three categories: ( 1 ) a link is 
critical if in its absence we must increase the number of 
driver nodes to maintain full control over the system. In 
this case the link is part of all maximum matchings of 
the network; ( 2 ) a link is redundant if it can be removed 
without affecting the current set of driver nodes (i.e. it 
does not appear in any maximum matching); (3) a link is 
ordinary if it is neither critical nor redundant (it appears 
in some but not all maximum matchings). Note that this 
classification can be efficiently done with a polynomial¬ 
time algorithm based on Berge’s property ( |Regin 1994), 
rather than enumerating all maximum matchings, which 
is infeasible for large networks. 


We can compute the density of critical (4 = Lc/L), 
redundant (Z,- = L^/L) and ordinary (Zq = Lq/L) links 
for a wide range of real-world networks. It turns out 
that most real networks have few or no critical links. 
Most links are ordinary, meaning that they play a role in 
some control configurations, but the network can be still 
controlled in their absence (Liu et al. 2011a). 


For model networks (ER and SF), we can calculate Zc, 
Zr, and Zo as functions of (k) (Fig. [T4|). The behavior of 
Zc is easy to understand: for small (k) all links are essen¬ 
tial for control (Zc « 1). As (fc) increases the network’s 
redundancy increases, decreasing Zc. The increasing re¬ 
dundancy suggests that the density of redundant links, 
Zf, should always increase with (fc), but it does not: it 
reaches a maximum at (Zc)c, after which it decays. This 
non-monotonic behavior results from a structural transi¬ 
tion driven by core percolation (Liu et al. 2011aI. Here, 
the core represents a compact cluster of nodes left in the 
network after applying a greedy leaf removal procedure: 
Recursively remove in-leaf (with Zcin = 1) and out-leaf 
(with fcout = 1 ) nodes’ neighbors’ all outgoing (or incom¬ 
ing) links. The core emerges through a percolation tran¬ 
sition (Fig. [^,d): for k < (Zc)c, Ucore = Ncore/N = 0, 
so the system consists of leaves only. At (fc)c a small 
core emerges, decreasing the number of leaves. For 
ER random networks, the analytical calculations predict 
(Zc)c = 2 e « 5.436564, in agreement with the numer¬ 
ical result (Fig. [T^), a value that coincides with (Zc) 
where Z^ reaches its maximum. Indeed, Z^ starts decaying 
at {k)c because after (Zc)^ the number of distinct max¬ 
imum matchings increases exponentially, which can be 
confirmed by calculating the ground state entropy using 
the cavity method (Liu et al.. 2011a). Consequently, the 
chance that a link does not participate in any control 
configurations decreases. For SF networks we observe 
the same behavior, with the caveat that (Zc)c decreases 
with 7 (Fig. [I4|:, d). 
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FIG. 14 (Color online) Link classification and core percola¬ 
tion. a, Dependence on (k) of the fraction of critical (red, 
Ic), rednndant (green, L) and ordinary (grey, G) links for an 
Erdos-Renyi (ER) network: peaks at (k) = {k}c = 2e and 
the derivative of is discontinuous at (fc) = (fe)c. b. Core 
percolation for the ER network occurs at (k) = (A:)c = 2e, 
which explains the h peak, c, d. Same as in a and b but for 
scale-free networks constructed using the static model. The 
ER and SF networks have N = 10^ nodes and the results are 
averaged over ten realizations with error bars defined as the 
standard error of the mean. Dotted lines are only a guide to 


the eye. After (Liu et al. 2011a|. 


2. Node Classification 



FIG. 15 (Color online) Emergence of bimodality in control¬ 
ling complex networks, (a) Ur and tic (insert) vs (k) in scale- 
free networks with degree exponents 7out = 7in = 3, display¬ 
ing the emergence of a bimodal behavior for high (fe). (b) 
Tir in scale-free networks with asymmetric in- and out-degree 
distribution, i.e. 7 c,ut = 3, 7in = 2.67 (upper branch) and 
7out = 2.67, 7in = 3 (lower branch). The control mode is 
pre-determined by their degree asymmetry. (c,d) Networks 
displaying centralized or distributed control. Both networks 
have Nu = 4 and A^c = 1 (red node), but they have rather 
different number of redundant nodes (blu e nodes). Nr = 23 
in (c) and Ar = 3 in (d). After (Jia et aL\ 20131. 


Given the existence of multiple driver node configura¬ 
tions, we can classify nodes based on their likelihood of 
being included in the minimum driver node set (MDNS): 
a node is (1) critical if that node must always be con¬ 
trolled to control the system, implying that it is part 
of all MDNSs; (2) redundant if it is never required for 
control, implying that it never participates in an MDNS; 
and (3) intermittent if it is a driver node in some control 


configurations, but not in others (Jia et al. 2013). 

For model networks with symmetric in- and out-degree 
distributions, we find that the fraction of redundant 
nodes (ur) undergoes a bifurcation at a critical mean de¬ 
gree {k)c'. for low (fc) the fraction of redundant nodes (rir) 
is uniquely determined by (k), but beyond (fc)c two dif¬ 
ferent solutions for coexist, one with very high and the 
other with very low value, leading to a bimodal behavior 


(Fig. 15 1 ). Hence for large (k) (after the bifurcation) two 
control modes coexist ( Jia et a^ [2013 ): (i) Centralized 
control: In networks that follow the upper branch of the 
bifurcation diagram most of the nodes are redundant, as 
in this case Ur is very high. This means that in these net¬ 
works only a small fraction of the nodes are involved in 
control (ric + ni is very low), hence control is guaranteed 
by a few nodes in the network. A good analogy would be 
a company involved in manufacturing whose leadership is 
concentrated in the hands of a few managers and the rest 
of the employees are only executors, (ii) Distributed con¬ 
trol: In networks on the lower branch Uc + n, can exceed 


90%. Hence, most nodes participate as driver nodes in 
some MDNSs, implying that one can engage most nodes 
in control. A good analogy would be an innovation-based 
horizontal organization, where any employee can take a 
leadership role, as the shifting tasks require. 

For ER random networks this bifurcation occurs at 
{k)c = 2e, corresponding to the core percolation thresh¬ 


old (Liu et al. 2012aI. 


Another way to assess a node’s importance for control 
is to quantify the impact of its removal on controllabil¬ 
ity. Consider a network with minimum number of driver 
nodes Ad. After a node is removed (deleted), denote 
the minimum number of driver nodes with A^. Once 
again, each node can belong to one of three categories: 
(1) A node is deletion critical if in its absence we have 
to control more driver nodes, i.e. A^ > Ad. For ex¬ 
ample, removing a node in the middle of a directed path 
will increase Ad. (2) A node is deletion redundant if 
in its absence we have A^ < Ad. For example, remov¬ 
ing a leaf node in a star will decrease Ad by 1. (3) A 
node is deletion ordinary if in its absence A^ = Ad. 
For example, removing the central hub in a star will not 
change Ad. The above node classification has been ap¬ 
plied to directed human protein-protein interaction net¬ 


works, whose directions indicate signal flow (Vinayagam 


et al.. 20161. In this context critical nodes tend to cor¬ 


respond to disease genes, viral tagets, through which a 
virus takes control over its host, and targets of FDA ap- 
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proved drugs, indicating that control-based classification 
can select biologically relevant proteins. 


3. Driver node classification 


To understand why a node is a driver node, we de¬ 


compose the driver nodes (IVd) into three groups (Ruths 


and Ruths, 2014): (1) source nodes (Ng) that have no 


incoming links, hence they must be directly controlled, 
being always driver nodes; (2) external dilations (iVe) 
arise due to a surplus of sink nodes (Nt) that have 
no outgoing links. Since each source node can con¬ 
trol one sink node, the number of external dilation is 
Ng = max(0, ~ Ns)] (3) internal dilations (Ni) oc¬ 
cur when a path must branch into two or more paths 
in order to reach all nodes (or equivalently a subgraph 
has more outgoing links than incoming links). This clas¬ 
sification leads to the control profile of a network de¬ 
fined as (?7s, ??e, ?7i) = {Ns/N,Ne/N,Ni/N), which quan¬ 
tifies the different proportions of control-inducing struc¬ 
tures present in a network. The measurements indicate 
that random network models do not reproduce the con¬ 
trol profiles of real-world networks and that the control 
profiles of real networks group into three well-defined 
clusters, dominated by external-dilations, sources, or 


internal-dilations (Ruths and Ruths, 2014). 


These results offer insight into the high-level organi¬ 
zation and function of complex networks. For exam¬ 
ple, neural and social network are source dominated, 
which allow relatively uncorrelated behavior across their 
agents and are thus suitable to distributed processing. 
Food webs and airport interconnectivity networks are 
internal-dilation dominated. They are mostly closed sys¬ 
tems and obey some type of conservation laws. In con¬ 
trast, trust hierarchies and transcriptional systems are 
external-dilation dominated. With their surplus sink 
nodes, these systems display correlated behavior across 
their agents that are downstream neighbors of a common 
source. 


G. Controllable Subspace, Control Centrality, and 
Structure Permeability 

Lin’s structural controllability theorem can tell us 
whether an LTI system (A,B) is structurally control¬ 
lable or not. If, however, the system is not structurally 
controllable, the theorem does not provide further infor¬ 
mation about controllability. Even if we are unable to 
make the system reach any point in the state space, we 
would like to understand which region of the state space 
is accessible to it, i.e., what region of the state space can 
we control it. For example, in the network of Fig. the 
control input ui is applied to the central hub Xi of the 
directed star with N = 3 nodes. The system is there- 



FIG. 16 (Color online) Control profiles of real and model net¬ 
works. The control profiles of real networks show a tendency 
to cluster around the three components {ris,r]e, rji) of the con¬ 
trol profile, implying that real networks broadly fall into three 
distinct classes: external-dilation dominated, source domi¬ 
nated, and internal-dilation dominated. The coloring of each 
small heatmap indicates the clustering observed in a wide 
range of real networks, with numbers in parentheses indicat¬ 
ing the number of networks present in each heatmap. Deeper 
shades of the heatmap represent a greater density of networks 


with control profiles located in that region. After (Ruths and 


Ruths 2014). 


fore stuck in the plane described by a 3 ia: 2 (t) = a 2 ia^ 3 (t), 
shaded in Fig.|^. Consequently, the network is not con¬ 
trollable in the whole state space, but it is controllable 
within the subspace defined by the plane. 


When we control a single node i, the input matrix B 
reduces to a vector b(j) with a single non-zero entry, and 
the controllability matrix C G becomes C{i). We 

can use rank(C(j)) as a natural measure of node i’s abil¬ 
ity to control the system. If rank(C(f)) = N, then node 
i alone can control the whole system. Any rank(C(i)) 
less than N yields the dimension of the subspace i can 
control. For example, if rank(C(z)) = I, then node i can 
only control itself. 


In reality the system parameters (i.e. the entries of 
A and B) are often not known precisely, except the ze¬ 
ros that mark the absence of connections, rendering the 
calculation of rank(C(j)) difficult. This difficulty can be 
again avoided using structural control theory. Assum¬ 
ing A and B are structured matrices, i.e., their elements 
are either fixed zeros or independent free parameters, 
then rank(C(i)) varies as a function of the free parame¬ 
ters of A and B. However, it achieves its maximum for 
almost all sets of values of the free parameters except 
for some pathological cases with Lebesgue measure zero. 


This maximal value is called the generic rank (Johnston 


et al.. 1984) of the controllability matrix C(i), denoted as 


rankg(C(z)), which also represents the generic dimension 
of the controllable subspace. 


We define the control capacity of a single node i, or 
control centrality, as the generic dimension of the con- 
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trollable subspace 


dc{A, b) = rankg C. 


(29) 


Here rankg C is the generic rank of the controllability 
matrix C associated with the structured system (A,b), 
where we control node i only (Liu et al. 2012bThis 


dehnition can also be extended to the case when we con¬ 
trol via a group of nodes. The above definition corre¬ 
sponds directly to our intuition of how powerful a single 
node is (or a group of nodes are) in controlling the whole 
network. For example, if the control capacity of a single 
node is iV, then we can control the whole system through 
it. 

The calculation of (ic(A, B) has a graph-theoretic in¬ 
terpretation (Hosoe 19801. Consider a structured sys¬ 


tem (A, B), in which all state vertices are accessible, and 
let us denote with Q the set of subgraphs of G(A, B) 
which can be spanned by a collection of vertex-disjoint 
cycles and stems. In this case, the generic dimension of 
the controllable subspace is 


dc(A,B) = max|F;(G)| 


(30) 


where |G(G)| is the number of edges in the subgraph G. 
This is called Hosoe’s controllable subspace theorem. Es¬ 
sentially, Hosoe’s theorem tells us that to calculate the 
generic dimension of the controllable subspace we need to 
find the cactus that contains as many edges as possible. 
Note that Hosoe’s theorem applies only to a structured 
system (A, B) that has no inaccessible state vertices. 
In calculating fic(A, B) for a general system (A, B), we 
should only consider the accessible part of the network. 

For a digraph with no directed cycles Hosoe’s theo¬ 
rem further simplifies: the controllability of any node 
equals its layer index: Cs{i) = h- Here the layer in¬ 
dex of a node is calculated from the unique hierarchi¬ 
cal structure of the digraph following a recursive labeling 


procedure (Liu et al. 2012b). For general networks, we 
can use linear programming to calculate dc(A,B) ( |PoI 
[iak[ 1990). We first get a new graph G'(A, B) from 


G(A,B) by adding to G(A,B) the edges {vi,V]s[+j) for 
i = 1, - ■ ■ ,N, j = 1, - ■ ■ ,M; and the loops {vi,Vi) for 
i = 1, • • • ,N + M, if they do not exist in G(A, B) (see 
Fig.[^. We associate the weight rce = 1 with every orig¬ 
inal edge e of G(A, B) and the weight We = 0 with every 
new edge. A collection of node-disjoint cycle in G'(A, B) 
covering all nodes will be called a cycle partition. It is 
easy to check that to calculate max^gg |G(G)| is equiv¬ 
alent to calculating the maximum-weight cycle partition 
in G'(A, B), which can then be solved by the follow¬ 
ing linear programming: rnax^ggg,,(-^ subject 

to: (1) ^ leaves node Vi) = 1 for every node 

Vi S G'(A, B); (2) '^{xg : e enters node Vi) = 1 for ev¬ 
ery node Vi S G'(A, B); (3) Xe € {0,1} for every edge 
eS G'(A,B). 



FIG. 17 (Color online) Calculation of control centrality, (a) 
The original controlled system is represented by a digraph 
G(A, B). (b) The modified digraph G'{A, B) used in solving 
the linear programming. Dotted and solid lines are assigned 
with weight Wij = 0 and 1, respectively. The maximum- 
weight cycle partition is shown in red, which has weight 3, 
corresponding to the generic dimension of controllable sub¬ 
space by controlling node x\ or equivalently the control cen¬ 
trality of node x\. 


Hosoe’s theorem also allows us to address a problem 
complementary to the notion of control centrality: iden¬ 
tify an optimal set of driver nodes of fixed cardinality 
M, denoted as Hd(AI), for a network of size N such that 
the dimension of the controllable subspace, denoted as 
\C{M)\, is maximized (Lo ludice et al. , 2015). If we solve 
this problem for each M G [1, A], we obtain a sequence 
of \C{M)\. To quantify the readiness or propensity of a 
network to be controllable, we can calculate the so-called 
network permeability measure (Lo ludice et al. 2015) 


/„^(|C(M)|-M)dM 


(31) 


Note that /r S [0,1]: 0 for N disconnected nodes, and 
1 for networks that are completely controllable by one 
driver node. Generally, for a network with a high per¬ 
meability, a large controllable subspace can be obtained 
with a reasonable small set of driver nodes. 


H. Controlling Edges 


So far we focused on nodal dynamics, where we moni¬ 
tored and controlled the state of nodes. The sole purpose 
of the edges was to pass information or influence between 
the nodes. In social or communication networks nodes 
constantly process the information received from their 
upstream neighbors and make decisions that are com¬ 
municated to their downstream neighbors. Most impor¬ 
tantly, in these systems nodes can communicate different 
information along different edges. Hence the information 
received and passed on by a node can be best represented 
by state variables defined on the incoming and outgoing 
edges, respectively. In this section we ask how to control 
systems characterized by such edge dynamics. 

To model s uch systems we place the state variables 
on the edges (Nepusz and Vicsek, 2012). Let y~(t) and 
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controllability of edge dynamics unchanged. 


(t) represent vectors consisting of the state variables ics (Liu et al. 2011a Posfai et al. 20131, but leaves the 
associated with the incoming and outgoing edges of node 
i, respectively. Let denote the fcout(*) x hnii) ma¬ 
trix. The equations governing the edge dynamics can be 
written as 


y+ (t) = M,y^ (t) - Ti 0 y+ (t) -h 


(32) 


where is a vector of damping terms associated with 
the outgoing edges, 0 denotes the entry-wise product of 
two vectors of the same size, and Ui = 1 if node z is a 
driver node and 0 otherwise. Note that even though the 
state variables and the control inputs are defined on the 
edges, we can still designate a node to be a driver node if 
its outgoing edges are directly controlled by the control 


I. Self-Dynamics and its Impact on Controllability 

The nodes of networked systems are often character¬ 
ized by some self-dynamics, e.g. a term of the form 
ii = aaXi, which captures the node’s behavior in the 
absence of interactions with other nodes. If we naively 
apply structural control theory to systems where each 
node has a self-dynamic term we obtain a surprising re¬ 
sult — a single control input can make an arbitrarily large 


inputs. Equation (32) states that the state variables of 2011a). This result represents a special case of the min- 


linear system controllable (Cowan et al ., 2012 Liu et al. 


the outgoing edges of node i are determined by the state 
variables of the incoming edges, modulated by a decay 
term. For a driver node, the state variables of its out¬ 
going edges will also be influenced by the control signals 
u^. Since each node i acts as a small switchboard-like 
device mapping the signals of the incoming edges to the 


outgoing edges using a linear operator M^, Eq. (32) is 
often called the switchboard dynamics. 

There is a mathematical duality between edge dynam¬ 
ics on a network G and nodal dynamics on its line graph 
£(G), which represents the adjacencies between edges of 
G. Each node of £{G) corresponds to an edge in G, and 
each edge in £(G) corresponds to a length-two directed 
path in G. By applying the minimum input theorem di¬ 
rectly to this line graph, we obtain the minimum number 
of edges we must drive to control the original network. 
However, this procedure does not minimize the number 
of driver nodes in the original network. This edge con¬ 
trol problem can be mapped to a graph theoretical prob¬ 


lem as follows (Nepusz and Vicsek, 2012). Define node 


i to be (i) divergent, if fcout(*) > fcin(*); (h) convergent, 
if kont{i) < kin{i); (iii) balanced, if kout{i) = kin{i). A 
connected component in a directed network is called a 
balanced component if it contains at least one edge and 
all the nodes are balanced. We can prove that the mini¬ 
mum set of driver nodes required to maintain structural 
controllability of the switchboard dynamics on a directed 
network G can be determined by selecting the divergent 
nodes of G and an arbitrary node from each balanced 
component. 

The controllability properties of this edge dynamics 
significantly differ from simple nodal dynamics. For ex¬ 
ample, driver nodes prefer hubs with large out-degree 
and heterogeneous networks are more controllable, i.e., 
require fewer driver nodes, than homogeneous net¬ 
works (Nepusz and Vicsek 2012). Moreover, positive 


correlations between the in- and out-degree of a node 
enhances the controllability of edge dynamics, without 


affecting the controllability of nodal dynamics (Posfai 


et al.: 2013). Conversely, adding self-loops to individ¬ 


imum input theorem: The self-dynamics contributes a 
self-loop to each node, hence each node can be matched 
by itself. Consequently, G(A) has a perfect matching, 
independent of the network topology, and one input sig¬ 


nal is sufficient to control the whole system (Liu et al. 


2011a). 


To understand the true impact of self-dynamics on net¬ 
work controllability, we must revisit the validity of the 
assumption that the system parameters are independent 
of each other. As we show next, relaxing this assumption 
offers a more realistic characterization of real systems, for 
which not all system parameters are independent. 

Assuming prototypical linear form of self-dynamics, 
e.g., first-order x = oqx, second-order x = gqX + aix, 
etc, we can incorporate the linear self-dynamics with the 
LTI dynamics of the network in a unified matrix form, 
as illustrated in Fig. |18| An immediate but counterin¬ 
tuitive result states that in the absence of self-dynamics 
riD is exactly the same as in the case when each node 
has a self-loop with identical weight w, i.e. each node is 
governed by precisely the same self-dynamics. This is a 
direct consequence of the identity 


rank[B,AB,--- ,A'''-^B] 
= rank[B, (A -f wljB, • • • , (A -|- wI)^-^B], 


(33) 


where on the left we have the rank of controllability ma¬ 
trix in the absence of self-loops, and on the right the same 
for a network where each node has an identical self-loop. 
For more general cases the minimum number of driver 


nodes Ad can be calculated from (27), i.e. the maximum 


geometric multiplicity of A’s eigenvalues. 

Note a remarkable symmetry in network controllabil¬ 
ity: If we exchange the fractions of any two types of self¬ 
loops with distinct weights, the system’s controllability. 


ual nodes enhances the controllability of nodal dynam- 


as measured by nu, remains the same (Fig. 19). For 
example, consider a network without self-loops. Equiv¬ 
alently, we can assume that each node contains a self¬ 
loop with weight zero. Then we systematically add more 
non-zero self-loops with identical weights to the network. 
Equivalently, we are replacing the zero-weight self-loops 
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Network topology Dynamic unit Integrated network Matrix form 



fln * 0 0 0 0 

0 «o 0 0 0 0 

0 0 flu * 0 * 

0 * 0 "o * 0 

0 * 0 0 Oq 0 

0 0 0 0 * a^' 



FIG. 18 (Color online) Integrating the network topology with 
nodal self-dynamics, (a) Ist-order self-dynamics: x = oqX. 
(b) 2nd-order self-dynamics: x = a^x -\- aix . (c) 3rd- 

order self-dynamics: x = aox + aix + a 2 X. To develop a 
graphical representation for the dth-order individual dynam¬ 
ics = aox^^^ + aix^^^ -I- • • • -I- ad-ix^^~^\ we denote each 
order by a colored square. The couplings among orders are 
characterized by links or self-loops. This graphical represen¬ 
tation allows the individual dynamics to be integrated with 
the network topology, giving rise to a unified matrix that re¬ 
flects the dynamics of the whole system. In particular, each 
dynamic unit in the unified matrix corresponds to a diagonal 
block and the nonzero elements (denoted by *) outside these 
blocks stand for the couplings among different dynamic units. 
Therefore, the original network of N nodes with order d self- 


dynamics is represented by a dN x dN matrix. After (Zhao 


et al. 20151. 


with non-zero self-loops, nu will first decrease as the frac¬ 
tion p of non-zero self-loops increases, reaching a mini¬ 
mum at p = After that, ud increases, reaching its 
maximum at p = 1, which coincides with n jy observed for 
p = 0 (Fig. 19 1 ). We can introduce more types of self¬ 
loops with different weights. If we exchange the fractions 
of any two types of self-loops, njy remains the same. This 
exchange-invariant property gives rise to a global symme¬ 
try point, where all the different types of self-loops have 
equal densities and the system displays the highest con¬ 
trollability (i.e., lowest number of driver nodes). This 
symmetry-induced optimal controllability holds for any 


network topology and various individual dynamics (Zhao 


et al. 2015). 


J. Control Energy 

Indentifying the minimum number of driver or actu¬ 
ator nodes sufficient for control is only the first step of 
the control problem. Once we have that, we need to ask 



FIG. 19 (Golor online) Impact of first-order individual dy¬ 
namics on the fraction of driver nodes njy. The values of the 
off-diagonal non-zero elements in A are randomly chosen and 
hence are independent, (a) nn in function of pa, the density 
of nodes that have the same type of nonzero self-loops. We 
observe a clear symmetry around ps = 1/2, indicating that 
ud reaches its minimum at pa = 1/2, where the densities of 
nodes with zero and non-zero self-loops are equal, (b) ud for 
an Erdos-Renyi random network with three types of self-loops 
Si, S2 and ss with densities pi^\ pP^ and pi^\ respectively. 
The color bar denotes the value of nu and the coordinates in 
the triangle stands for pi^^ pP^ and pP\ There is a global 
symmetry point where the three types of self-loops have the 
same density 1/3, and no reaches its minimum value. After 
(Zhao et al. 2015 I. 


an equally important question: How much effort is re¬ 
quired to control a system from a given set of nodes? 
The meaning of the term “control effort” depends upon 
the particular application (Kirk 2004a| ). In the case of a 
rocket being thrust upward, the control input u{t) is the 
thrust of the engine, whose magnitude \u(t)\ is assumed 
to be proportional to the rate of fuel consumption. In or¬ 
der to minimize the total expenditure of fuel, the control 
effort can be defined as |M(<)|dt, which is related to 
the energy consumed by the rocket. In the case of a volt¬ 
age source driving a circuit containing no energy storage 
elements, the source voltage is the control input u(t) and 
the source current is directly proportional to u(t). If the 
circuit is to be controlled with minimum energy dissipa¬ 
tion, we can define the control effort as fg u^(t)dt, which 
is proportional to the energy dissipation. If there are sev¬ 
eral control inputs, the general form of control effort can 
be defined as u'^(t)R(t)u(t) dt, where R(t) is a real 
symmetric positive-definite weighting matrix. 


Consider the LTI system (11) driven from an arbitrary 


initial state Xi towards a desired final state Xf by the 
external signal u(t) in the time interval t G [0,r]. We 
define the associated control effort in the quadratic form 

£-(T)^ r||u(t)f dt, (34) 

Jo 


called the “control energy” in the literature (Chen et al. 


2015 Yan et al. 2012 2015). Note that (34) may not 
have the physical dimension of energy, i.e., M T“^, 
in real control problems. But for physical and electronic 
systems we can always assume there is an hidden con- 
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stant in the right-hand side of (34) with proper dimen¬ 


sion, which ensures that £{T) has the dimension of en¬ 
ergy. In many systems, like biological or social systems. 


where (34) does not correspond to energy, it captures the 


effort needed to control a system. 

For a fixed set of driver nodes the control input u(t) 
that can drive the system from Xi to Xf can be chosen 
in many different ways, resulting in different trajectories 
followed by the system. Each of these trajectories has its 
own control energy. Of all the possible inputs, the one 
that yields the minimum control energy is 

u{t) = exp(A'^(T - t))W-i(T)vf, (35) 
where W (t) is the gramian matrix 


W(t) = / exp(AT)BB'^ exp(A"'’r) dr, (36) 


which is nonsingular for any i > 0 (Lewis et al. 


2012). Note that W(oo) is known as the controllabil¬ 
ity Gramian, often denoted with Wc ( Kailath[ 1980). 


The energy associated with the optimal input (35) is 


£{T) = v/W~^(T)vf, where Vf = Xf — exp(Ar)xi rep¬ 
resents the difference between the desired state under 
control and the final state during free evolution without 
control. Without loss of generality, we can set the final 
state at the origin, Xf = 0, and write the control energy 
as 


£(r) = x/H-^(T)xi 


(37) 


where H(T) = exp(—AT)W(T) exp(—A"'"T) is the sym¬ 
metric Gramian matrix. We can further define the nor¬ 
malized control energy as 


E{T) = 


£{T) 


x^H 


(38) 


When Xi is parallel to the direction of one of H’s eigen¬ 
vectors, the inverse of the corresponding eigenvalue corre¬ 
sponds to normalized energy associated with controlling 
the system along the particular eigen-direction. 

Using the Rayleigh-Ritz theorem, the normalized con¬ 
trol energy obeys the bounds 


Vmax — ^min 


< E{T) < 


(39) 


where and the maximum and minimum eigen¬ 


values of H, respectively (Yan et al. 2012). 


Assuming linear individual dynamics characterized by 
the self-loop an = — (a -I- Si) where Si = the 

strength of node i and a is a parameter that can make the 
symmetric A (describing an undirected network) either 
positive or negative definite, we can choose a single node 
with index c as the driver node. In this case, the lower 


and upper energy bounds follow 

'T-i 
E, 


[(A+AT)-i], 

T -1 ^ 0 


small T 

large T, A is PD 
large T, A is semi PD 


,exp {2Xjs[T) —0 large T, A is not PD 


(40) 


Er, 


'T-® ( 0 > 1 ) 

= e(A, c) 

T-i ^ 0 

^exp (2AiT) —0 


Here Ai > At > 


> 


small T 

large T, A is not ND 
large T, A is semi ND 
large T, A is ND 

(41) 

Aat are the eigenvalues of A, 


and e(A, c) is a positive energy that depends on the ma¬ 
trix A and the choice of the controlled node c. PD (or 
ND) means positive-definite (or negative-definite), re¬ 
spectively. The scaling laws (40) and (41) can be gen¬ 


eralized to directed networks, in which case the decay 
exponents Ai and Xn are replaced by ReAi and RcAat, 
respectively. 


Equations (40) and (411 suggest that the scaling of the 


control energy is rather sensitive to the control time T. 
For small T, in which case we wish to get our system very 
fast to its destination, both Emin and Umax decay with 
increasing T, implying that setting a somewhat longer 
control time requires less energy. For large T, however, 
we reach a point where we cannot reduce the energy by 
waiting for longer time. This occurs when the system has 
its equilibrium point in the origin, then any attempt to 
steer the system away from the origin must overcome a 
certain energy barrier. 

The control energy is rather sensitive to the direction of 


et al. 


the state space in which we wish to move the system (Yan 
2015|). To see this, consider a scale-free network 


with degree exponent 7 . If we drive the system through 
all its nodes (Ad = A), the control energy spectrum, 
describing the probability that moving in a randomly 
chosen eigen-direction will require energy £, follows the 
power law P{£) ^ £~'^. Consequently, the maximum 
energy required for control depends sublinearly on the 
system size, fmax ~ implying that even in the 

most costly direction the required energy grows slower 
than the system size. In other words, if we control each 
node, there are no significant energetic barriers for con¬ 
trol. If, however, we aim to control the system through 
a single node (Ad = 1 ), the control spectrum follows a 
power law with exponent —1, i.e., P{£) ~ £~^, which 
only weakly depends on the network structure. There¬ 
fore the maximum energy required for control increases 
as fmax e^. This exponential increase means that 
steering the network in some directions is energetically 
prohibitive. Finally, if we drive a finite fraction of nodes 
(1 < Ad < A), the control spectrum has multiple peaks 
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x(() = AK{t) + Bu(t) 



0 0.5 1 1.5 2 2.5 3 

t 



FIG. 20 (Color online) Energy spectrum, (a) A three-node 
weighted network can be controlled via a single control input 
u(t), injected to the driver node shown in red. The input 
matrix B is reduced to a vector (1,0,0)"*". Each node has a 
negative self-loop, which makes all eigenvalues of the state 
matrix A negative, hence stable, (b) The optimal control 
signals that minimize the energies required to steer the net¬ 
work from the initial state xq = x(0) = (0,0, 0)^ to three 
different desired states Xd = x(t) at t = 3, with the con¬ 
straint ||xd|| = 1. (c) The trajectories of the network state 
x(t) driven by the control inputs shown in (b). (d) The energy 
surface for all normalized desired states, i.e., |jxd|| = 1, which 
is an ellipsoid spanned by the controllability Gramian’s three 
eigen-directions (arrows). The ellipsoid nature of the spec¬ 
trum illustrates the widely different energies we need to move 
the network shown in (a) in different directions in the state 
space. The squares correspond to the three cases depicted in 


(b) and (c). After (Yan et al. 20151. 


and the maximum energy required for control scales as 
i^max ~ . Hence, as we increase the number of 

driver nodes, the maximum energy decays exponentially. 

These results raise an important question: in case 
of 1 < Ad < N, how to choose the optimal set of 
Ad driver nodes such that the control energy is min¬ 
imized? Such a combinatorial optimization problem 
(also known as the actuator placement problem) has 
not been extensively studied in the literature. Only re¬ 
cently has it been shown that several objective func¬ 
tions, i.e. energy-related controllability metrics associ¬ 
ated with the controllability Gramian Wc of LTI sys¬ 
tems (e.g. Tr(W“^), log(detWc), rank(Wc)), are actu¬ 
ally submodular (Cortesi et al.\ 2014[ Summers et al. 


2014|^Summers and Lygeros[ 2013). A submodular func¬ 


tion [j / has the so-called diminishing returns property 
that the difference in the function value that a single el¬ 


ement X makes when added to an input set X decreases 
as the size of the input set increases. The submodularity 
of objective functions allows for either an efficient global 
optimization or a simple greedy approximation algorithm 
with certain performance guarantee to solve the combi¬ 
natorial optimization problems (Nemhauser et al. 1978). 
In particular, the submodularity of those energy-related 
controllability metrics has been explored to address the 
actuator placement problem in a model of the European 


power grid ( 

Cortesi et al. 

2014 

Summers et al. 

2014 

Summers and Lygeros 201? 

)■ 


K. Control Trajectories 


So far we have focused on the minimization of 
driver/actuator nodes and the energy cost of controlling 
LTI systems. The characteristics of the resulting con¬ 
trol trajec tories are also in t eresti ng and worthy of ex¬ 
ploration (Sun and Motter 2013). A state x*^°) of the 


LTI system is called strictly locally controllable (SLC) if 
for a ball centered at x^^^ with radius £ > 0 

there is a constant d > 0 such that any final state 
inside the ball can be reached from x(°) with a 

control trajectory entirely inside the ball B{x^^\e) (see 
Fig. 21 1 ). Figure [2T|3 shows that in a two-dimensional 
LTI system xi = xi + ui{t), ±2 = xi, for any state in 
the Xi > 0 half-plane, the minimal-energy control tra¬ 
jectories to any neighboring final state with a smaller 
a: 2 -component will necessarily cross into the Si < 0 half¬ 
plane. 

It has been shown that for a general LTI system when¬ 
ever the number of control inputs is smaller than the 
number of state variables (i.e., Ad < A), then almost all 
the states are not SLC (Sun and Motter [2013 ). There¬ 
fore, the minimal-energy control trajectory is generally 
nonlocal and remains finite even when the final state is 
brought arbitrarily close to the initial state. The length 
|lx(t)||dt of such a trajectory generally increases with 
the condition number of the Gramian. Furthermore, the 


optimal control input (351 that minimizes the energy 


cost ||u(t)|pdt will fail in practice if the controlla¬ 
bility Gramian (36) is ill conditioned. This can occur 


even when the controllability matrix is well conditioned. 
There is a sharp transition, called the controllability tran¬ 
sition, as a function of the number of control inputs, be¬ 
low which numerical control always fails and above which 
it succeeds. These results indicate that even for the sim¬ 
plest LTI dynamics, the disparity between theory and 
practice poses a fundamental limit on our ability to con¬ 


trol large networks (Sun and Motter, 2013). 


® Denote V{S) as the power set (i.e. the set of all the subsets) of a K that satisfies f{X U {a;}) — f{X) > f{y U {a:}) — /(T), for any 

set S. Then a submodular function is a set function / :E{cS)— >■ X C y C S and x e cS \ T. 
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FIG. 21 (Color online) Strictly local controllability, (a) Illus¬ 
tration of a state that is strictly locally controllable (left) and 
a state that is not (right), (b) The state space of a simple 
LTI system with two state variables xi and X 2 - The curves 
indicate control trajectories of minimal energy for given ini¬ 
tial state (open symbol) and final states (solid symbols). The 
arrows in the background represent the vector field in the ab¬ 
sence of control input ui (t). Note that any state that is not on 
the line a^i = 0 is not SLC, because the minimal-energy con¬ 
trol trajectories to any neighboring final state with a smaller 
2 : 2 -component will necessarily cross into the xi < 0 half-plane. 
After | |Sun and Motter] |2013| ). 


Indeed, we usually don’t use the minimum-energy con¬ 
trol input (35) to steer the system to desired final states, 
simply because it is an open-loop (or non-feedback) con¬ 
troller 1^ which tends to be very sensitive to noise. A 
more practical and robust strategy is to use a simple 
linear feedback control to bring the system asymptoti¬ 
cally towards a certain state, while minimizing the energy 
cost. This is a typical objective of optimal control the¬ 
ory, which aims to design control signals that will cause 
a process to satisfy some physical constraints and max¬ 
imize (or minimize) a chosen performance criterion (or 
cost function) ( Kirk[ 2004b Naidu 2002). 


III. CONTROLLABILITY OF NONLINEAR SYSTEMS 

So far we focused on the controllability of linear sys¬ 
tems. Yet, the dynamics of most real complex systems 
is nonlinear, prompting us to review the classical results 
on nonlinear controllability and their applications to net¬ 
worked systems. 

Consider a control system of the form 

x = f(x,u), (42) 

where the state vector x is in a smooth connected man¬ 
ifold JH of dimension N, and the control input u € U 


^ An open-loop control system does not use feedback. The control 
input to the system is determined using only the current state 
of the system and a model of the system, and is totally indepen¬ 
dent of the system’s output. In contrast, in a closed-loop control 
system, the output has an effect on the input (through feedback) 
so that the input will adjust itself based on the output. 


is a subset of Note that (42) has been frequently 
used to model the behavior of physical, biological and 


social systems (Hermann and Krener, 1977). Roughly 


speaking, (42) is controllable if one can steer it from any 


point Xq S Ad to any other point Xi G Ad by choosing u 
from a set of admissible controls U, which is a subset of 
functions mapping M+ to lA. 

The controllability of nonlinear systems has been ex¬ 


tensively studied since the early 1970s (Brocket! 1972 


Conte et al. 200^ Ell^ 1970 de Figueiredo and Chen 


1993 

Haynes and Hermesl 

1970[ IHermann and Krener 

1977 

Isidori 

1995[ Lobry[ 

1970 INijmeijer and van der 

SchaftI 1990 

Rugh[ 1981| ISontag 1998 Sussmann and 


Jurdjevic 1972). The goal was to derive results of similar 


reach and generality as obtained for linear time-invariant 
systems. However, this goal turned out to be too am¬ 
bitious, suggesting that a general theory on nonlinear 
controllability may not be feasible. Fortunately, as we 
discuss in this section, the concerted effort on nonlinear 
control has led to various weaker notions of nonlinear 
controllability, which are easier to characterize and often 
offer simple algebraic tests to explore the controllability 
of nonlinear systems. 


A. Accessibility and Controllability 


As we will see in the coming sections, we can rarely 
prove or test controllability of an arbitrary nonlinear sys¬ 
tem. Instead, we prove and test weaker versions of con¬ 
trollability called local accessibility and local strong ac¬ 
cessibility. We start by defining these notions. 

Accessibility concerns the possibility to reach or access 
an open set of states in the state space from a given ini¬ 


tial state. If the system (42) is locally accessible from an 
initial state Xq then we can reach or access the neigh¬ 
borhood of Xq through trajectories that are within the 


neighborhood of Xq. Mathematically, the system (42) 


is called locally accessible from Xq if for any non-empty 
neighborhoods V C Ad of Xg and any ti > 0, the reach¬ 
able set Ti^{xQ,< ti) contains a non-empty open set. 
The system is called locally accessible if this holds for 
any Xg. Here, the reachable set 77.^(xg, < ti) includes all 
states that can be reached from Xg within a time ti, fol¬ 
lowing trajectories that are within the neighborhood of 
Xg. Mathematically, the reachable set from xg in time up 
to ti is defined as 7^^(xg, < ti) = U,-<tj7^^(xg, r). Here 
77 .^(xo,t) is the reachable set from Xg at time r > 0 
following trajectories that remain in V for t < t. 

If we look at states that can be reached exactly at time 
H, then we have a stronger version of local accessibility. 


System (42) is said to be locally strongly accessible from 


Xg if at any small time H > 0 the system can reach or 
access the neighborhood of Xg through trajectories that 
are within the neighborhood of xg. Mathematically, this 
means that for any non-empty neighborhoods V of Xg and 
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any ti > 0 sufficiently small, the reachable set 7^^(xo, ti) 
contains a non-empty open set. If this holds for any 
Xq, then the system is called locally strongly accessible. 
Clearly, local strong accessibility from xg implies local 
accessibility from Xg. The converse is generally not true. 

Local controllability asks whether the system is con¬ 
trollable in some neighborhood of a given state. Math¬ 
ematically, the system (42) is called locally controllable 
from Xg if for any neighborhood V of Xg, the reachable 
set 7?.^(xg, < ti) is also a neighborhood of Xg for any ti 
small enough. The system is called locally controllable if 
this holds for any Xg. Clearly, local controllability im¬ 
plies local accessibility. It turns out that for a large class 
of systems local controllability implies local strong acces¬ 
sibility. But the converse is not always true. 

If we do not require the trajectories of the system to re¬ 
main close to the starting point, i.e., we allow excursions, 
then we have the notion of global controllability. System 
(42) is globally controllable from Xg if the reachable set 
from Xg is M itself, i.e., 7^(xg) = Ut^>g7^^(xg, ti) = M. 
In other words, for any xi € A4, there exists ti > 0 and 
u : [0,ti] — ?> U such that the solution of ( |4^ starting at 
Xg at time 0 with control u{t) satishes x(ti) = Xi. If this 
holds for all Xg G AI, then the system is called globally 
controllable. 

Complete algebraic characterizations of global control¬ 
lability of nonlinear systems have proved elusive. Weaker 
notions of controllability are easier to characterize than 
controllability. For example, it can be proven that for 
some nonlinear systems, accessibility can be decided in 


polynomial time, while controllability is NP-hard (Son- 


tag 1988). For complex networked systems we expect 


that only weaker notions of controllability can be char¬ 
acterized. 


If the linearized control system is controllable (in the 
sense of a linear time-invariant system), then for any 
e > 0 the original nonlinear system is locally control¬ 
lable from X*, where the control functions u(-) are taken 
from the set Ue. 

In other words, many real systems operate near some 
equilibrium points and in the vicinity of such points, con¬ 
trollability can be decided using the tools developed for 
linear systems, discussed in the previous sections. 


2. Linearization around a trajectory 


We can also study the linearized control system along 
a trajectory. Consider a nonlinear control system in the 


form of (42). A trajectory represents the path the system 


follows as a function of time in the state space. It can be 
mathematically defined as a function (x, u) : [Tg,ri] —>■ 

X 


O, where O is a nonempty open subset of 
and x(t 2 ) = x(ti) -|- f(x(t), u(t)) dt, for all £ 


[7g,Ti]. The linearized control system of (42) along a 
trajectory (x, u) : [rg,Ti] —>• O is a linear time-varying 
control system x = A(t)x -|- B(t)u with t £ [Tq, Ti], and 


Mt) = ^(x(0:U(t)), B{t) = |^(x(t),u(t)). (44) 


If the linearized control system along the trajectory 
(x,u) : [rg,Ti] —)■ O is controllable in the sense of a 
linear time-varying system, then the original nonlinear 
system is locally controllable along the trajectory. Once 
again, this means that we can use linear control theory 
to explore the controllability of nonlinear systems. 


B. Controllability of Linearized Control System 


It is typically difficult to test the controllability of a 
nonlinear system. Yet, as we discuss next, studying the 
controllability properties of its linearization around an 
equilibrium point or along a trajectory can often offer 


an efficient test of local nonlinear controllability (Coron 


2009). 


1. Linearization around an equilibrium point 


Consider an equilibrium point (x*, u*) G M xU of the 
nonlinear control system (42), meaning that f(x*,u*) = 
0. Assume that U contains a neighborhood of u*. For 
e > 0, we define a set of control functions Ue = {u(-) G 
U||lu(t) -u*|| < e,t > 0}. The linearized control system 
at (x*, u*) is a linear control system x = Ax -f Bu with 

di di 

A=^(x*,u*), B=-(x*,u*). (43) 


3. Limitations of linearization 


The linearization approaches described above may 
sound powerful, but they have severe limitations. First, 
they only provide information about controllability in the 
immediate vicinity of an equilibrium point or a trajec¬ 
tory. Second and most important, it may be the case 
that the linearized control system is not controllable, but 
the original nonlinear system is actually controllable. 


Consider, for example, a model of a front-wheel drive 
car with four state variables: the positions (a;i,a; 2 ) of 
the center of the front axle, the orientation (j) of the 
car, and the angle 9 of the front wheels relative to the 
car orientation (Fig. [2^. There are two control inputs 
{ui,U 2 ), where ui, the steering velocity, represents the 
velocity with which the steering wheel is turning, and U 2 
is the driving velocity. Assuming that the front and rear 
wheels do not slip and that the distance between them is 


I = 1, the car’s equations of motion have the form (Nel- 
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FIG. 22 (Color online) Controlling a car. The figure shows 
a model of a front-wheel drive car with 4 state variables 
{x\,X2,4>y&) and 2 control inputs (iti,U 2 ). While this system 
is globally controllable (see Sec. III.El, its linearized dynam- 
ics around the origin is not controllable. Figure redrawn from 


(Sontag 19981. 


son 


1967 Sontag 1998) 


(il\ 


h 


^cos{9 + (f)'^ 

X2 

= Ui 

0 

0 

+ U2 

sin(0 -1- (f) 

(/> 



sin 6 

[^j 


vJ 


\ 0 / 


(45) 


will focus on control-affine systems, referring the reader 
to (Hermann and Krener 1977 Sontag, 1998) for more 
general nonlinear systems. 


C. Basic concepts in differential geometry 

Before we discuss the nonlinear tests for accessibility 
and controllability, we need a few concepts in differential 
geometry, like Lie brackets and distributions. 


1. Lie brackets 


For nonlinear control systems, both controllability and 
accessibility are intimately tied to Lie brackets. The rea¬ 
son is simple. In the nonlinear framework, the directions 
in which the state may be moved around an initial state 
Xq are those belonging to the Lie algebra generated by 
vector fields f(xo,u), when u varies in the set of admis¬ 
sible controls U ( Isidori[ 1995 Sontag 1998). Here the 
Lie algebra A generated by a family F of vector fields is 
the set of Lie brackets [f.g] with f, g S and all vec¬ 
tor fields that can be obtained by iteratively computing 
Lie brackets. In turn, a Lie bracket is the derivative of a 
vector field with respect to another. 

Consider two vector fields f and g on an open set D C 
The Lie bracket operation generates a new vector 
field [f,g], defined as 


The linearization of (45) around the origin is 




'^M2^ 

±2 


0 



0 

vv 


\Ul) 


(46) 


which is uncontrollable, because X 2 and (j) are time- 
invariant and not controlled by any of the system’s in¬ 
puts. Yet, from our driving experience we know that 
a car is controllable. We will prove that this system is 
indeed globally controllable in Sec. |HLE[ 

System (45) belongs to an especially interesting class 
of nonlinear systems, called control-affine systems, where 
f(x, u) is linear in the control signal u 


M 

X = f(x) 4-^gi(x)M*. (47) 

2 = 1 


[f, g](x) = 

where ^ and ^ are the Jacobian matrices of g and f, 
respectively. Higher order Lie brackets can be recursively 
defined as 


ad^g(x) = g(x), (49) 

adfg(x) = [f,adf"^g](x), Vfc > 1 (50) 

where “ad” denotes “adjoint”. 

To understand the physical meaning of the Lie bracket, 
consider the following piece-wise constant control inputs 


u(t) 


01,0)'^, te[o,T) 

1(0,1)'^, tG[r,2r) 

I tG[2r,3T) 

[(0,-1)'^, tG[3r,4T) 


(51) 


Here, f is called the drift vector field, or simply drift', 
and gi, • • • , gM are called the control vector fields. The 


system (47) is called driftless if f(x) = 0, which arises in 
kinematic models of many mechanical systems, e.g., in 


(45). Control-affine systems are natural generalization 


of linear time-invariant systems. Many nonlinear con¬ 
trollability results were obtained for them. Hereafter, we 


applied onto a two-inputs control-affine system 

X = gi(x)ai-kg2(x)M2 


(52) 


with initial state x(0) = Xq ( Sastry[ 1999). The piece- 
wise constant control inputs (51) can be considered as a 


sequence of “actions” applied for example to a car (gi. 
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FIG. 23 Lie bracket. The physical meaning of a Lie bracket 
can be demonstrated by applying the piece-wise constant con¬ 
trol inputs ( |51[ ) to a two-inputs driftless control-affine system 
X = gi(x)ui -I- g2(x)u2. Up to terms of order r^, the differ¬ 
ence between the final state x( 4 r) and the initial state xo is 
given by the Lie bracket [gi,g2](xo). 


g 2 , reverse-gi, reverse-g 2 ). In the limit r 
state reached at t = 4r is 


0 the final 


x(4r) = Xo -h ^^gi(xo) - ^g 2 (xo) j -b 0{t^). 

(53) 

We see that up to terms of order the state change is 
exactly along the direction of the Lie bracket [gi, g2](xo) 


dgi 


(see Fig. 23). 


Consider two examples that demonstrate the meaning 
of Lie brackets. First, the Brockett system is one of the 
simplest driftless control-afRne systems (iBrockett 1982) 


Xi{t) = Ui 
X2{t) = U2 
isit) = U2X1 - U1X2 


(54) 


which can be written in the form of (52) using gi (x) = 
(1, 0, —a;2)"'’ and g2(x) = (0,1, Xi)"^, or equivalently gi = 
(^-X 2 j^,andg 2 = These two operators 

gi and g2 have a nontrivial Lie bracket [gi,g2](x) = 
g 3 (x) = 2(0,0,1)"”", or equivalently ga = Consider 


the system ( [54| initially at the origin, hence gi( 0 ) = 
(1, 0,0)"'", g2(0) = (0,1, 0)"^. If we again apply the control 
sequence ( |51[ ) with time interval r = 1, we can check 
that the final state reached at t = 4 is (0,0,2)"'", which is 
precisely captured by [gi,g2](0). 

Note that for the Brockett system we have 
[gi, [gi,g2]](x) = [g2, [gi,g2]](x) = 0 . A similar three- 
dimensional Lie algebra, called the Heisenberg algebra, 
also arises in quantum mechanics. Hence the Brockett 


system is also known as the Heisenberg system (Bloch 


2003). Note, however, that the commutation relations 


obeyed by the Heisenberg algebra do not always ap¬ 
ply to general nonlinear systems. To see this consider 
again the model of a front-wheel drive car ( [4^ , repre¬ 
senting another two-input control-affine system, where 
the two control vector fields gi = (0,0,0,1)"”" and g 2 = 
(cos(0-|-(^),sin(0-|-(()),sin0,O)'^ can be interpreted as the 


actions steer and drive, respectively. Some Lie brackets 
from gi(x) and g2(x) are 


g3(x) = [gi,g2](x) = 


g 4 (x) = [[gl,g 2 ],g 2 ](x) = 


^ — sin(0 -b (/))^ 

cos(6> -b (j)) 


cos 9 

V 

0 J 

^ — sin 

'l — 

coscj) 

) — 

0 


1 0 / 


(55) 


(56) 


Equation (55) can be interpreted as [steer, drive] = 
wriggle, arising from the sequence of actions 
(steer, drive, reverse steer, reverse drive), which is 
what we do in order to get a car out of a tight 


parking space. Similarly, (56) can be interpreted as 


[wriggle, drive] = slide, arising from the sequence of 
actions (wriggle, drive, reverse wriggle, reverse drive), 
which is what we do during parallel parking. Equations 


(55 p6| indicate that starting from only two control 


inputs: steer and drive, we can “generate” other actions, 
e.g., wriggle and slide, which allows us to fully control 
the car. 

The above two examples demonstrate that by apply¬ 
ing the right sequence of control inputs we can steer the 
system along a direction that the system does not have di¬ 
rect control over. In general, by choosing more elaborate 
sequences of control inputs we can steer a control-afEne 
system in directions precisely captured by higher-order 
Lie brackets, e.g., [g2, [gi,g2]], [[gi,g2], [g2, [gi,g2]]], etc. 
If the system of interest has a drift term f, we also have 
to consider Lie brackets involving f. This is the reason 
why nonlinear controllability is closely related to the Lie 
brackets. 


2. Distributions 

To discuss the nonlinear tests of accessibility and con¬ 
trollability, we need the notion of distribution in the sense 
of differential geometry. A distribution can be roughly 
considered as the nonlinear version of the controllability 
matrix of a linear system. 

Consider m vector fields gi, g2, • ’ ’ ; Sm on an open set 
V c K^. We denote 

A(x) = span{gi(x),g 2 (x), • • • ,gr„(x)} (57) 

as the vector space spanned by the vectors 
gi(x),g 2 (x),-• • ,g,„(x) at any fixed x e X>. Es¬ 
sentially, we assign a vector space A(x) to each point 
X in the set V. The collection of vector spaces A(x), 
X S I? is called a distribution and referred to by 

A = span{gi,g 2 , • • • ,gm}- 


(58) 
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If the vectors gi(x),g 2 (x), • • • ,gm(x) are linearly inde¬ 
pendent for any x in V, then the dimension of A(x) is 
constant and equals m. In this case we call A a non¬ 
singular distribution on V. For example, in the Brockett 
system we have gi(x) = (1,0,g 2 (x) = (0,1, 
g 3 (x) = [gi,g 2 ](x) = (0,0,2)'^. Since gi(x),g 2 (x),g 3 (x) 
are linearly independent for all x G we conclude 
that the distribution A = spanjgi, g 2 , gs} is nonsingu¬ 
lar. Similarly, in the front-wheel drive car system of 
Fig. gi(x),g 2 (x),g 3 (x) and g 4 (x) are linearly in¬ 
dependent for all X S hence the distribution A = 
spanjgi,g 2 ,g 3 ,g 4 } is nonsingular. Note that a nonsin¬ 
gular distribution is analogous to a full rank matrix. 


{f,gi,-- - ,gM} and iteratively generating new, inde¬ 
pendent vector fields using Lie brackets. This can be 
achieved by constructing the Philip Hall basis of the Lie 
algebra, which essentially follows a breadth-first search 
and the search depth is defined to be the number of 


nested levels of bracket operations (Duleba 1998 Serre 


1992). 


In general, accessibility does not imply controllability, 
which is why accessibility is a weaker version of control¬ 
lability. Consider a simple dynamical system 



(62) 


D. Nonlinear Tests for Accessibility 

1. Accessibility 


Roughly speaking, accessibility concerns whether we 
can access all directions of the state space from any given 
state. The accessibility of control-affine systems can be 
checked using a simple algebraic test based on Lie brack¬ 
ets. 

For control-affine systems (471, we denote C as the lin¬ 
ear combinations of recursive Lie brackets of the form 


[Xfe,[Xfc_i,[... ,[X2 ,Xi]...]]],/c = 1,2,... , (59) 

where X^ is a vector field in the set {f, gi, • • • , gu}- As 
the linear space C is a Lie algebra, it is closed under 
the Lie bracket operation. In other words, [f,g] e C 
whenever f and g are in C. Hence C is called as the 
accessibility algebra. 

The accessibility distribution C is the distribution gen¬ 
erated by the accessibility algebra C: 

C'(x) = span{X(x)|XeC}. (60) 

Consider a control-affine system ^T\ and a state Xq G 
M C If 


dimC'(xo) = N 


(61) 


then the system is locally accessible from Xq. Equation 


(61) is often called the aecessibility rank condition (ARC) 
at Xg. If it holds for any Xg, then the system is called 
locally accessible. 

Interestingly, the sufficient ARC is “almost” necessary 
for accessibility. Indeed, if the system is accessible then 
AR C holds for a l l x in an op en and dense subset of 
(|lsidori| |1995] [Sontagi [iMl. 


The computation of the accessibility distribution C 
is nontrivial, because it is not known a priori how 
many (nested) Lie brackets of the vector fields need 
to be computed until the ARC holds. In practice, a 
systematic search must be performed by starting with 


which can be written in the control-affine form (47) with 
f(x) = (a:|,0)'^ and g(x) = (0,1)"'’. We can compute 
some Lie brackets: [f,g](x) = —(2x2,0)'^, [f, [f, g]](x) = 
(2,0)’’’. Since [f, [f,g]](x) is independent from g(x), we 
conclude that dimC'(x) = 2, for any state in K^, indicat¬ 
ing that the system is locally accessible. But the system 
is not locally controllable: > 0 for all X 2 ^ 0, 

i.e., xi always grows as long as the system is not at the 
a; 2 -axis. In other words, the drift vector field f always 
steers the system to the right unless X 2 = 0. 

If we compute the accessibility distribution C for a lin¬ 
ear system x = Ax -f Bu = Ax -|- where B = 

[bi, • • • ,bM], we find that C'(xg) is spanned by Axg to¬ 
gether with the constant vector fields b^, Ab^, A^b^, • • •, 
for * = 1, • • • , M. More precisely, 

C'(xo) = span{Axg} + Im(B, AB, A^B, • • • , A^-^B), 

(63) 

where Im() stands for the image or column space of a 
matrix. Note that the term span{Axg} does not appear 


in Kalman’s controllability matrix (12). Only at Xg = 0, 


Eq.(63) reduces to Kalman’s controllability matrix. This 


shows that accessibility is indeed weaker than controlla¬ 
bility, because the former does not imply the latter while 
the latter induces the former. 


2. Strong accessibility 

A nonlinear test for strong accessibility tells us whether 
we can reach states in the neighborhood of the ini¬ 
tial state exactly at a given small time. Define Cg as 
the strong accessibility algebra, i.e., the smallest algebra 
which contains gi,g2,’-- , gM and satisfies [f, w] e Cg, 
Vw S Cg. Note that Cg C C and Cg does not contain 
the drift vector field f. Define the corresponding strong 
accessibility distribution 

C'g(x) = span{X(x)|X S Cg}. (64) 

If dimC'g(xg) = N then the system is locally strongly 
accessible from Xg. If this holds for any Xg, then the 
system is called locally strongly accessible. If we compute 
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the strong accessibility distribution C for a linear system 
(A,B), we will find that 


C'o(xo) = Im(B, AB, A^B, • • • , A^-^B). (65) 

Then dimC'o(xo) = A is equivalent with Kalman’s rank 
condition ( |T^ . In other words, strong accessibility and 
controllability are equivalent notions for linear systems. 


E. Nonlinear Tests for Controllability 


For general nonlinear systems, we lack conditions that 
are both sufficient and necessary for controllability. Yet, 
as we discuss next, we have some sufficient conditions 
that are believed to be almost necessary as well. 

Consider a special class of control-affine system ( [T^ 
with f(x) e span{gi(x), • • • ,gM{x)} for all x G A1 C 
In other words, the drift vector field f(x), which 
describes the intrinsic dynamics of the system, can be 
spanned by the control vector fields gi(x),--- ,gM{x)- 
Then, if dimC'(xo) = A, the system is locally controllable 
from Xq. If this holds for all x G Af, then the system is 
globally controllable. 

Driftless systems (f(x) = 0 ), like the front-wheel drive 
car system (45), naturally fall into this class. To see 


this, we recognize that the determinant of the matrix 
formed by the vectors gi(x), g 2 (x), g 3 (x) = [gi,g 2 ](x) 
and g 4 (x) = [[gi,g 2 ],g 2 ](x), i.e.. 


det 


(q cos[9 + (j)) — sin(0-I-())) —sincf^ 
0 sin(0 -|- (j)) cos{9 + cj)) cos (j) 

0 sin 9 cos 9 0 

yi 0 0 0 y 


( 66 ) 


is identically equal to 1, regardless of x, implying that 
dimC(xo) = A = 4 for all Xq G Hence the front- 
wheel drive car systems is globally controllable, in line 
with our physical intuition and experience. 

For control-affine systems that do not fall into the 
above two classes, Sussmann provided a general set of 
sufficient conditions (Sussmann 1987). We call a Lie 


bracket computed from {f, gi,--- ,gM} bad if it con¬ 
tains an odd number of f factors and an even number 
of each g^ factors. Otherwise we call it good. The de¬ 
gree of a bracket is the total number of vector fields from 
which it is compuated. Denote with permuta¬ 
tion group on M symbols. For a G ^ ^ 

bracket computed from {f,gi,--- ,gM}, define tt(b) as 
the bracket obtained by fixing f and changing g^ by gcr(fe), 
1 < k < M. The control-afhne system (47) is locally con¬ 


trollable from Xq if dimC(xo) = A and every bad bracket 
b has the property that ,0(b)(xo) = X^o-gSm is 

a linear combination of good brackets, evaluated at Xq, 
of degree lower than b. 




Motif 6 


Motif 8 


icSp) 


FIG. 24 (Color online) Symmetries and controllability. The 
eight different three-node neuronal network motifs studied in 
( [Whalen et al.\ |2015| ) . Those motifs display a variety of sym¬ 
metries. For example. Motif 1 has a full S 3 symmetry, and 
Motif 3 has a reflection S 2 symmetry across the plane through 
node 2. Not all symmetries have the same effect on network 
controllability. 


F. Controllability of Nonlinear Networked Systems 


1 . Neuronal network motifs 


While most complex systems are described by nonlin¬ 
ear continuous-time dynamics defined over a network, 
there has been little attention paid so far to the controlla¬ 
bility of such systems, due to obvious mathematical chal¬ 
lenges. Controllability studies of continuous-time non¬ 
linear dynamics are still limited to very simple networks 
consisting of a few nodes, like neuronal network motifs 


governed by Fitzhugh-Nagumo dynamics (Whalen et al. 


2015). These offered an opportunity to study the im¬ 


pact of structural symmetries on nonlinear controllabil¬ 
ity. The three-node neuronal motifs shown in Fig.p4|can 
have multiple symmetries. Yet, not all symmetries have 
the same effect on network controllability. For example, 
with identical nodal and coupling parameters. Motif I has 
a full S3 symmetry, rendering the poorest controllability 
over the entire range of coupling strengths. Similarly, no 
controllability is obtained from node 2 in Motif 3, which 
has a reflection S 2 symmetry across the plane through 
node 2. Surprisingly, the rotational C3 symmetry in Mo¬ 
tif 7 does not cause loss of controllability at all. Note 
that symmetries have an impact on network controllabil¬ 
ity in linear systems as well. For example, in the case of a 
directed star with LTI dynamics for which we control the 
central hub (Fig. |^, a symmetry among the leaf nodes 
renders the system uncontrollable. 

Extending this analysis to larger networks with sym¬ 
metries remains a challenge, however. Group represen¬ 
tation theory might offer tools to gain insights into the 
impact of symmetries on the controllability of nonlinear 


networked systems (Whalen et al. 2015). Note, however. 


that for large real networks such symmetries are less fre¬ 
quent. 
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2. Boolean networks 


The controllability of Boolean networks, a class of 
discrete-time nonlinear systems that are often used 
to model gene regulations, has been intensively stud¬ 


ied (Akutsu et al. 2007 Cheng and Qi 20091. We can 


prove that hnding a control strategy leading to the de¬ 
sired final state is NP-hard for a general Boolean network 
and this problem can be solved in polynomial time only 
if the network has a tree structure or contains at most 


one directed cycle (Akutsu et al. 20071. Interestingly, 


based on semi-tensor product of matrices (Cheng et al. 


20071 and the matrix expression of Boolean logic, the 


Boolean dynamics can be exactly mapped into the stan¬ 


dard discrete-time linear dynamics (Cheng and Qi, 20091. 


Necessary and sufficient conditions to assure controllabil¬ 


ity of Boolean networks can then be proved (Cheng and 


Qi 20091. Despite the formally simplicity, the price we 


need to pay is that the size of the discrete-time linear dy¬ 
namical system is 2^, where N is the number of nodes in 
the original Boolean network. Hence, the controllability 
test will be computationally intractable for large Boolean 
networks. 


IV. OBSERVABILITY 

Before controlling a system, it is useful to know its po¬ 
sition in the state-space, allowing us to decide in which 
direction we should steer it to accomplish the control ob¬ 
jective. The position of a system in the state-space can be 
identified only if we can measure the state of all compo¬ 
nents separately, like the concentration of each metabo¬ 
lite in a cell, or the current on each transmission line of a 
power grid. Such detailed measurements are often infea¬ 
sible and impractical. Instead, in practice we must rely 
on a subset of well-selected accessible variables (outputs) 
which can be used to observe the system, i.e. to estimate 
the state of the system. A system is said to be observable 
if it is possible to recover the state of the whole system 
from the measured variables inputs and outputs). This 
is a fundamental and primary issue in most complex sys¬ 
tems. 

In general, we can observe a system because its com¬ 
ponents form a network, hence the state of the nodes 
depend on the state of their neighbors’. This offers the 
possibility to estimate all unmeasured variables from the 
measured ones. If the inputs and model of the system are 
known, observability can be equivalently defined as the 
possibility to recover the initial state x(0) of the system 
from the output variables. 

To be specihc, let us assume that we have no knowledge 
of a system’s initial state x(0), but we can monitor some 
of its outputs y{t) in some time interval. The observabil¬ 
ity problem aims to establish a relationship between the 
outputs y{t), the state vector x(t), and the inputs u(t) 


such that the system’s initial state x(0) can be inferred. 
If no such relation exists, the system’s initial state can¬ 
not be estimated from the experimental measurements, 
i.e., the system is not observable. In other words, if the 
current value of at least one state variable cannot be de¬ 
termined through the outputs sensors, then it remains 
unknown to the controller. This may disable feedback 
control, which requires reliable real-time estimates of the 
system’s state. 

Note that observability and controllability are math¬ 
ematically dual concepts. Both concepts were first in¬ 
troduced by Rudolf Kalman for linear dynamical sys¬ 


tems (Kalman 1963), and were extensively explored in 


nonlinear dynamical systems by many authors ([Besanon 


2007 Diop and Fliess 1991a|b[ Hermann and Krener 


1977 Isidori 1995 Sontag and Wang, 1991). 


In this section, we first discuss methods that test the 
observability of linear and nonlinear control systems. We 
also discuss the parameter identifiability problem, which 
is a special case of the observability problem. Finally, 
we introduce a graphical approach to identify the min¬ 
imum set of sensor nodes that assure the observability 


of nonlinear systems (Aguirre and Letellier 


iKhan 


and Doostmohammadian |2011| 

Khan and Moura 

12008 

Letellier and Aguirre 2005| Lete 

her et al.\ 2006 Letellier 

and Aguirre 

2010 

Siddhartha and van Schuppen 

2001) 


and its application to metabolic networks (Liu et al. 


2013). 


A. Observability Tests 

1. Linear systems 


For linear systems there is an exact duality between 
controllability and observability. To see this, consider an 
LTI control system 


x(t) = Ax(t) -t- Bu(t) 
y{t) = Cx(t). 


(67a) 

(67b) 


The duality principle states that an LTI system (A, B, C) 
is observable if and only if its dual system (A"'’, B"^) 

is controllable. Mathematically, the duality can be 
seen and proved from the structure of the controllability 
Gramian and the observability Gramian. In terms of net¬ 
work language the duality principle has a straightforward 
interpretation: The linear observability of a network A 
can be addressed by studying the controllability of the 
transposed network A"”", which is obtained by flipping 


the direction of each link in A (Fig. 25). 


Thanks to the duality principle, many observability 
tests can be mapped into controllability tests. For exam¬ 
ple, according to Kalman’s rank condition, the system 
(A,B,C) is observable if and only if the observability 
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b 


X, 

X 2 

X 3 


FIG. 25 (Color online) Duality principle. If a system follows 
the LTI dynamics (67a I, the observability of the network A 



shown in (a) can be addressed by studying the controllability 
of the transposed network shown in (b), obtained by re¬ 
versing the direction of each link. This is a general property 
of all networks. 


matrix 


O = 


c 

CA 

CA2 


CA 


N-1 


( 68 ) 


has full rank, i.e., rankO = N (Kalman, 1963 Luen- 


berger} 19791. This rank condition is based on the fact 


that if the N rows of O are linearly independent, then 
each of the N state variables can be determined by linear 
combinations of the output variables y(t). 


2. Nonlinear systems 

Consider a nonlinear control system with inputs u(t) € 
and outputs y{t) G 


x(t) = f(t,x(t),u(t)) 
y{t) = h{t,x{t),u{t)) 


(69) 


where f(-) and h(-) are some nonlinear functions. 

Mathematically, we can quantify observability from ei¬ 


ther an algebraic viewpoint (Conte et al ., 2007 Diop 


and Fliess 1991a|b 

) or a d 

point ( 

Hermann and Krener 


19771. Here we focus on 


the former. If a system is algebraically observable, then 
there are algebraic relations between the state variables 
and the successive derivatives of the system’s inputs and 
outputs (Diop and Fliess 1991a|b ). These algebraic rela¬ 
tions guarantee that the system is observable and forbid 
symmetries. A family of symmetries is equivalent to in¬ 
finitely many trajectories of the state variables that fit 
the same specified input-output behavior, in which case 
the system is not observable. If the number of such tra¬ 
jectories is finite, the system is called locally observable. 
If there is a unique trajectory, the system is globally ob¬ 
servable. 


Consider, for example, the dynamical system defined 
by the equations 


Xi = X 2 X 4 -j- u 

±2 = X2X3 
±3=0 
X4 = 0 

y = xi 


(70) 


The system has a family of symmetries ay. 
{xi, 0 : 2 , ^ 3 , 0 : 4 } —?> {xi,\x 2 ,X 3 ,X 4 /X\, so that the 
input u and the output y and all their derivatives 
are independent of A ( |Anguelova[ 2004 1 . This means 
that we cannot distinguish whether the system is in 
state (cci,a; 2 , 3 : 3 , 3 : 4 )"^ or its symmetric counterpart 
{xi,Xx 2 ,X 3 ,X 4 /X)'^, because they are both consistent 
with the same input-output behavior. Hence we cannot 
uncover the system’s internal state by monitoring xi 
only. 

The algebraic observability of a rational system is de¬ 
termined by the dimension of the space spanned by the 
gradients of the Lie-derivatives 


^ ~ dx, 

i=l 


Ell". 

jgn 1=1 


( 1 + 1 ) 


d 


du\ 


U) 


(71) 


of its output functions h(t, x(t), u(<)). The observability 
problem can be further reduced to the so-called rank test: 


the system (69) is algebraically observable if and only if 


the NM X N Jacobian matrix 


J = 


dL°jh 

dxi 


dL)hM 

dx\ 


dx\ 


dV;-'- 

dx\ 


dL°f>n 

dx2 


dL°hM 

dx2 


d^hi 

dxN 


dL°hM 


dx]\ 


dv: 


dL^: 


dx2 


dx]s 






dx2 


dx]\ 


(72) 


has full rank (Diop and Fliess 1991a|b|, i.e.. 


rank J = N. 


(73) 


Note that for an LTI system ( 67a|67b ), the Jacobian ma¬ 
trix (72) reduces to the observability matrix (68). 


For rational dynamic systems, the algebraic observabil¬ 
ity test can be performed using an algorithm developed 


by Sedoglavic (Sedoglavic 2002). The algorithm offers a 


generic rank computation of the Jacobian matrix (72) us¬ 


ing the techniques of symbolic calculation, allowing us to 
test local algebraic observability for rational systems in 
polynomial time. This algorithm certifies that a system 
is locally observable, but its answer for a non-observable 
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system is probabilistic with high probability of success. 
A system that is found non-observable can be further 
analyzed to identify a family of symmetries, which can 
confirm the system is truly non-observable. 


B. Minimum sensor problem 


In complex systems, the state variables are rarely in¬ 
dependent of each other. The interactions between the 
system’s components induce intricate interdependencies 
among them. 

Hence a well-selected subset of state variables can con¬ 
tain sufficient information about the remaining variables 
to reconstruct the system’s complete internal state, mak¬ 
ing the system observable ( |Liu et aL\ |2013| ). 

We assume that we can monitor a selected subset of 
state variables, i.e. y(t) = (• • ■ , • • •)"'’, correspond¬ 

ing to the states of several nodes that we call sensor 
nodes or just sensors. Network observability can then 
be posed as follows: Identify the minimum set of sen¬ 
sors from whose measurement we can infer all other state 
variables. For linear systems, this problem can be solved 
using the duality principle and solving the minimum in¬ 
put problem of the transposed network For general 


nonlinear systems this trick does not work. While (73) of 


fers a formal answer to the observability issue and can be 
applied to small engineered systems, it has notable prac¬ 
tical limitations for large and complex systems. First, it 
can only confirm if a specific set of sensors can be used 
to observe a system or not, without telling us how to 
identify them. Therefore, a brute-force search for a min¬ 


imum sensor set requires us to inspect via (731 about 
2^ sensor combinations, a computationally prohibitive 
task for large systems. Second, the rank test of the Jaco¬ 
bian matrix via symbolic computation is computationally 


limited to small systems (Sedoglavic 2002). Hence, the 


fundamental question of identifying the minimum set of 
sensors through which we can observe a large complex 
system remains an outstanding challenge. 

To resolve these limitations, we can exploit the 
dynamic interdependence of the system’s components 


through a graphical representation ( 

Khan and Doostmo- 

hammadian 20111 Linl 19741 Murotal 2009 

Reinschke 

1988 Siddhartha and van Schuppenj 2001). The proce- 

dure consists of the following steps 

Liu et al. 

2013 

I: 


(i) Inference diagram: Draw a directed link Xi —)■ xj 
if Xj appears in xfs differential equation (i.e., if is 
not identically zero), implying that one can retrieve some 
information on Xj by monitoring Xi as a function of time. 
Since the constructed network captures the information 
flow to infer the state of individual variables, we call it 
the inference diagram (Fig. [26j:). 

(ii) Strongly connected component (SCC) decomposi¬ 
tion: Decompose the inference diagram into a unique set 


of maximal SCCs (dashed circles in Fig. 26:), i.e. the 


3 Reactions 


A + B + C —► D + F + J 
D-<-^E 
H + I 

J + K —►G + H 


Xi 

i2 

X 4 

X5 

X6 

X7 

Xs 

X 9 

^10 

Xll 



Balance equations 

= —kiXiX2X3 

— ^k^x^x^x^ 

— 3723/3 

= +^1X1X2373 - k 2 X 4 + ^3X5 

= +k2X4 - ksX5 

= +^ 1 X 1 X 2273 

= +^4X3X9 - ^5X7 + fc 6 a 7 ioa 7 ii 

= -kiXsXQ + ^5X7 + kexioxn 
~ ^4X3X9 “1“ /35X7 

= +fciXiX 2 X 3 — fesXiQXii 
= —fcsXioXii 





non-root SCC 

6 

root SCC 

• 

sensor node 

A 0 

non-sensor nod^ 


FIG. 26 (Color online) The graphical approach to determine 
the minimum sensors of a chemical reaction system, (a) A 
chemical reaction system with eleven species (A,B,- • ■, J,K) 
involved in four reactions. Since two reactions are reversible, 
we have six elementary reactions, (b) The balance equations 
of the chemical reaction system shown in (a). The concentra¬ 
tions of the eleven species are denoted by xi,X2, ■ ■ ■ , ®io, xn, 
respectively. The rate constants of the six elementary reac¬ 
tions are given by k\,k 2 ,--- ,ke, respectively. The balance 
equations are derived using the mass action kinetics, (c) The 
inference diagram is constructed by drawing a directed link 
{xi —> Xj) as long as Xj appears in the RHS of Xi’s bal¬ 
ance equation shown in (b). Strongly connected components 
(SCCs) are marked with dashed circle. Root SCCs, which 
have no incoming links, are shaded in grey. A potential min¬ 
imum set of sensor nodes, whose measurements allows us to 
reconstruct the state of all other variables (metabolite con¬ 
centrations), are shown in red. After (Liu et aT] 2013). 


largest subgraphs chosen such that in each of them there 
is a directed path from every node to every other node 
(Cormen et al. 1990). Consequently, each node in an 
SCC contains some information about all other nodes 
within the SCC. 

(iii) Sensor node selection: Those SCCs that have no 
incoming edges are referred to as root SCCs (shaded cir¬ 
cles in Fig. 26:). We must choose at least one node from 
each root SCC to ensure the observability of the whole 
system. For example, the inference diagram of Fig. |26[ : 
contains three root SCCs; hence we need at least three 
sensors to observe the system. 

The graphical approach (GA) described above can be 
used to determine whether a variable provides full observ- 
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ability of small dynamic systems (Aguirre et al. 2008 
Letellier and Aguirre, 2005| ) . As these systems have only 
a few state variables, steps (ii) and (iii) are often not 
necessary. For large networked systems, the GA is very 
powerful because it reduces the observability issue, a dy¬ 
namical problem of a nonlinear system with many un¬ 
knowns, to a property of the static graph of the infer¬ 
ence diagram, which can be accurately mapped for an 
increasing number of complex systems, from biochemical 
reactions to ecological systems. 

We can prove that monitoring the root SCCs identi¬ 
fied by the GA are necessary for observing any nonlinear 


dynamic system (Liu et al. 20131. In other words, the 


number of root SGGs yields a strict lower bound for the 
size of the minimum sensor set. Gonsequently, any state 
observer (i.e. a dynamical device that aims to estimate 
the system’s internal state) will fail if it doesn’t monitor 
these sensors. 

If the dynamics is linear, the duality principle maps 
the minimum sensor problem into the minimum input 
problem and predicts not only the necessary, but also 
the sufficient sensor set for observability. Numerical sim¬ 
ulations on model networks suggest that for linear sys¬ 
tems the sufficient sensor set is noticeably larger than the 


necessary sensor set predicted by GA (Liu et al. 2013). 
This is because that any symmetries in the state variables 
leaving the inputs, outputs, and all their derivatives in¬ 


variant will make the system unobservable (Sedoglavic 


2002 1 . For structured linear systems, the symmetries cor¬ 


respond to a particular topological feature, i.e., dilations, 
which can be detected from the inference diagram. Yet, 
for general nonlinear systems, the symmetries can not be 
easily detected from the inference diagram only. 

For linear systems the minimum sensor set predicted 
by the GA is generally not sufficient for full observabil¬ 
ity. Yet, for large nonlinear dynamical systems the sym¬ 
metries in state variables are extremely rare, especially 
when the number of state variables is big, hence the sen¬ 
sor set predicted by GA is often not only necessary but 


also sufficient for observability (Liu et al. 2013). 


To better understand network observability, next we 
apply the developed tools to biochemical and technolog¬ 
ical networks. 


1. Biocfnemical reaction systems 

Consider a biochemical reaction system of N 
species { 5 i, 52 ,--- ^Sn} involved in R reactions 
{ 7 ^ 1 , ^ 2 , * • ■ with 


N 


N 


where aji > 0 and j3ji > 0 are the stoichiometry coeffi¬ 


cients. For example, (74) captures the reaction 2 H 2 


O 2 = 2 H 2 O with ail = 2 , ai 2 = 1 and Pn = 2 . 


Under the continuum hypothesis and the well-mixed 


assumption the system’s dynamics is described by (69), 


where Xi{t) is the concentration of species Si at time 
t, the input vector u(t) represents regulatory signals 
or external nutrient concentrations, and the vector y(t) 
captures the set of experimentally measurable species 
concentrations or reaction fluxes. The vector v(x) = 
(ni(x), i; 2 (x), • • • ,'!;ii;(x))'^ is often called the flux vec¬ 


tor, which follows the mass-action kinetics (Heinrich and 
Schuster[ 1996[ Palsson, 2006) 


N 


= %!!■ 


(75) 


with rate constants kj > 0. The system’s dynamics is 
therefore described by the balance equations 


ii = /*(x) = ^riji;j(x), 

i=i 


(76) 


where Fi, = Pji — aji are the elements of the N x R 


stoichiometric matrix F. The RHS of (76) represents a 


sum of all fluxes Vj that produce and consume the species 

5 ,. 

Assuming that the outputs y{t) are just the concen¬ 
trations of a particular set of sensor species that can 
be experimentally measured, then observability prob¬ 
lem aims to identify a minimum set of sensor species 
from whose measured concentrations we can determine 
all other species’ concentrations. In this context, the ad¬ 
vantage of GA is that it does not require the system’s 
kinetic constants (which are largely unknown in vivo), re¬ 
lying only on the topology of the inference diagram. For 
a metabolic network or an arbitrary biochemical reaction 
system, the topology of the inference diagram is uniquely 
determined by the full reaction list, which is relatively 


accurately known for several model organisms (Schellen- 


berger et al. 2010). Applying GA to biochemical reaction 


systems offers several interesting results, elucidating the 


principles behind biochemical network observability (Liu 


et al.. 2013): 


a) Species that are not reactants in any reaction, being 
instead pure products, will be root SCCs of size one. Con¬ 
sequently, they are always sensors, and must be observed 


by the external observer (e.g., xq in Fig. 26:). 


b) For root SCCs of size larger than one (e.g. { 0 : 4 , 0 : 5 } 
and {xr, xs, xg} in Fig. 26:), any node could be chosen as 
a sensor. Given that some root SCCs are quite large, and 
typically we only need to monitor one node for each root 
see, the number of sensor nodes is thus considerably 
reduced. 

c) A minimum set of sensors consists of all pure prod¬ 
ucts and one node from each root SCC of size larger than 


one (e.g. {x;i„xg,xr} in Fig. 26:). 


d) Since any node in a root SCC can be selected as 
a sensor node, there are Hg = equivalent 
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sensor node combinations, representing the product of all 
root SCCs’ sizes. For example, in Fig. [26] : we have three 
root SCCs with sizes rii = 1,2,3, hence fig = 1x2x3 = 6. 
This multiplicity offers significant flexibility in selecting 
experimentally accessible sensors. 


It turns out that the minimum set of sensors obtained 
by GA almost always achieve full observability for the 


whole system, except in some pathological cases (Liu 


et al. 20131. The sufficiency of the sensors predicted 


by GA is unexpected because substantial details about 
the system’s dynamics are ignored in GA, hence offering 
an exact proof that the sufficiency of the predicted sen¬ 
sors for observability is a difficult, if not an impossible, 
task. Note, however, that the rigorous proof of sufficiency 
and the systematic search for exceptional cases making a 
system unobservable remain open questions. 


2. Power grid 


In the power grid, the state variables represent the 
voltage of all nodes, which in practice can be determined 
by phasor measurement units (PMUs). Since a PMU 
can measure the real time voltage and line currents of 
the corresponding node, a PMU placed on a node i will 
determine the state variables of both node i and all of 
its first nearest neighbors. In this case the observabil¬ 
ity problem can be mapped to a purely graph theoretical 
problem. The random placement of PMUs leads to a net¬ 
work observability transition (Yang et al. 2012), which 
is a new type of percolation transition that characterizes 
the emergence of macroscopic observable components in 
the network as the number of randomly placed PMUs in¬ 
creases (Fig. 27). Using the generating function formal¬ 
ism (Newman et al. 2001), we can analytically calculate 
the expected size of the largest observable component for 
networks with any prescribed degree distribution. This 


has been demonstrated for real power grids (Yang et al. 


2012). Moreover, it has been found that the percolation 


threshold decreases with the increasing average degree or 
degree heterogeneity (Yang et al. 2012). 


The random placement of PMUs apparently will not 
solve the minimum sensor problem. For a power grid, the 
problem of identifying the minimum set of sensor nodes is 
reduced to the minimum dominating set (MDS) problem: 
Identify a minimum node set D C V for a graph G = 
{V,E) such that every node not in D is adjacent to at 
least one node in D (Fig. [28^,b). Gonsider a undirected 
network G. Node i is either empty (with occupation state 
Ci = 0) or occupied by sensors (with a = 1). In other 
words, if Ci = 1 then node i can be considered a sensor 
node. Node i is called observed if it is a sensor node 
itself or it is not a sensor node but adjacent to one or 
more sensor nodes. Otherwise node i is unobserved. The 
MDS problem requires us to occupy a minimum set D of 



Fraction of directly observed nodes, </> 


FIG. 27 (Color online) Observability transitions in the power 
grid, (a) Fraction of the largest observable component as a 
function of the fraction of directly observed nodes (</>) in net¬ 
works with prescribed degree distributions of the power grids 
of Eastern North America (black), Germany (red), Europe 
(green), and Spain (blue). The continuous lines are analytical 
predictions, and the symbols represent the average over ten 
10®-node random networks for 10 independent random PMU 
placements each. The inset shows a magnification around 
the transitions, with the analytically predicted thresholds (j>c 
indicated by arrows. After (Yang et al. 20121. 


nodes so that all N nodes of G are observed. 0 

The MDS problem for a general graph is NP-hard, and 
the best polynomial algorithms can only offer dominat¬ 
ing sets with sizes not exceeding log N times of the mini¬ 


mum size of the dominating sets (Lund and Yannakakis 


1994 Raz and Safra, 1997). If the underlying network 


has no core, we can solve exactly the MDS problem in 
polynomial time using a generalized leaf-removal (GLR) 
process (Fig.[2^,d). The GLR process can be recursively 
applied to simplify the network G. If eventually all the 
nodes are removed, then the set of nodes occupied during 
this process must be an MDS and choosing them as sen¬ 


sor nodes will make the whole network observable (Zhao 


et al. 2014). If, however, the final simplified network is 


non-empty, then there must be some nodes that are still 
unobserved after the GLR process. The subnetwork in¬ 
duced by these unobserved nodes is referred to as the core 
of the original network G. For networks with an exten¬ 
sive core, a belief-propagation algorithm, rooted in spin 
glass theory, can offer nearly optimal solutions, which 


also performs well on real-world networks (Zhao et al. 


2014). Recently, probabilistic methods have been devel¬ 


oped to approximate the size of the MDS in scale-free 
networks (Molnar et al. 2014). 


® Interestingly, the MDS problem can also be formalized as a con¬ 
trol problem on a undirected network by assuming that every 
edge in a network is bi-directional and every node in the MDS 
can control all of its outgoing links separately | |Jose and Tatsuy^ 
|2012||. This formulation has recently been applied to analyze bi¬ 
ological networks |Nacher and Akutsu[|201^ Wuchty[|2014[|. 
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FIG. 28 (Color online) Dominating set and generalized leaf 
removal process, (a-b) Dominating set. A dominating set 
of a graph G = {V, E) is a subset D oi V such that every 
vertex not in D is adjacent to at least one vertex in D. A 
minimum dominating set (MDS, shown in blue) is a domi¬ 
nating set of the smallest size, (c-d) Generalized leaf removal 
(GLR) process. If a network is sufficiently sparse, then its 
minimum dominating set (MDS) can be found exactly using 
GLR, consisting of two basic operations illustrated in (c) and 
(d). Blue circles denote nodes occupied with sensor nodes. 
White circles denote empty (i.e. non-occupied) and unob¬ 
servable nodes. Green circles denote empty but observable 
nodes, (c) For an empty leaf node i, its only adjacent node 
j must be occupied, i.e. be chosen as a sensor node. Conse¬ 
quently all adjacent nodes of j are observed. Node j and its 
adjacent nodes can be removed from the network to simplify 
the MDS problem, (d) If an empty observed node i has only 
a single unobserved adjacent node j, then it must be an opti¬ 
mal strategy not to occupy node i. Hence, the link between i 
and j can be removed from the network to simplify the MDS 


problem. After (Zhao et al. 20141. 


to Xt in the inference diagram. For example, in Fig. |26| ;, 
X 4 can only be inferred from X 5 while xi can be inferred 
from any other nodes, b) There are important differ¬ 
ences in the complexity of the inference process, which 
depends on the size of the subsystem we need to infer 
for a given sensor choice. The SCO decomposition of the 
inference diagram indicates that to observe Xt from Xs, 
we need to reconstruct A 4 = gS metabolite con¬ 
centrations, where 5s denotes the set of all SCCs that 
are reachable from Xg, and is the size of the i-th SCO. 
This formula can be extended to multiple targets, c) To 
identify the optimal sensor node for any target node, we 
can minimize c 5 which is the minimum amount 
of information required for the inference process. For ex¬ 
ample, if Xt is inside an SCC of size larger than one (e.g., 
xi in Fig. 26:), then the optimal sensor can be any other 
node in the same SCC (e.g., X 2 or X 3 in Fig. 26:). If all 
other nodes in the same SCC is experimentally inaccessi¬ 
ble, then the optimal sensor node belongs to the smallest 
SCC that points to Xi (e.g., Xg in Fig. [26]:). Note that 
this minimization procedure can be implemented for any 
inference diagram in polynomial time. Hence the graph¬ 
ical approach can aid the efficient selection of optimal 
sensors for any targeted node, offering a potentially in¬ 
dispensable tool for biomarker design. 


D. Observer Design 


The observability test and the graphical approach men¬ 
tioned above do not tell us how to reconstruct the state 
of the system from measurements. To achieve this we 
must design an observer, a dynamic device that runs a 
replica of the real system, adjusting its state from the 
available outputs to uncover the missing variables. 


For an LTI system (67a 67b), we can easily design the 


so-called Luenberger observer (Luenberger, 1964 1966 


1971) 


C. Target observability 

In many applications it is overkill to observe the full 
system, but it is sufficient to infer the state of a subset 
of target variables. Such target variables could for ex¬ 
ample correspond to the concentrations of metabolites 
whose activities are altered by a disease (]Barabasi et al. 


2011), representing potential biomarkers. In case those 
target variables cannot be directly measured, we can in¬ 
voke target observability, and aim to identify the optimal 
sensor(s) that can infer the state of the target variables. 
These could represent the optimal experimentally acces¬ 
sible biomarkers for a disease. The graphical approach 
discussed above helps us select such optimal sensors: a) 
The state of a target node Xt can be observed from a 
sensor node Xs only if there is a directed path from Xg 


z{t) = A z(t) -f L [y(t) — C z(t)] -)- B u(t) (77) 

where the N x K matrix L is to be specihed later. Note 
that with initial condition z(0) = x(0), the Luenberger 
observer will follow z{t) = x(t) exactly for all t > 0. 
Because x(0) is typically unaccessible, we start from 
z(0) ^ x(0) and hope that z{t) will asymptotically con¬ 
verge to x(t), i.e. the state of the observer tracks the state 
of the original system. This can be achieved by choos¬ 
ing a proper L matrix such that the matrix [A — L C] 
is asymptotically stable, in which case the error vec¬ 
tor e(t) = z(t) — x(t), satisfying e(t) = [A —LC]e(t), 
will converge to zero with rate determined by the largest 
eigenvalue of [ A — L C]. 

For nonlinear systems the observer design is rather in¬ 


volved and still an open challenge (Besangon 2007 Fried- 


land 1996). 
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1. Parameter Identification 


Most modeling efforts assume that the system param¬ 
eters, like the rate constants of biochemical reactions, 
are known. Yet, for most complex systems, especially in 
biological context, the system parameters are usually un¬ 
known or are only known approximately. Furthermore, 
the known parameters are typically estimated in vitro, 
and their in vivo relevance is often questionable. This 
raises a natural question: Can we determine the model 
parameters through appropriate input/output measure¬ 
ments, like monitoring the concentrations of properly se¬ 
lected chemical species? This problem is called param¬ 


eter identification (PI) in control theory (Bellman and 


Astroml 119701 |Glad and Ljungl 

1990 

ILjung 1987 Poh- 

janpalol 1978 Saccomani et al.\ 

2003 

1 . 


We can formalize the parameter identifiability prob¬ 
lem as the observability problem of an extended sys¬ 
tem as follows ( Anguelova[ 2004). For this we consider 
the system parameters 0 as special state variables with 
time-derivative zero (d0/dt = 0). We can extend the 
state vector to include a larger set of state variables, i.e., 
(x(t), 0), allowing us to formally determine whether/how 
the system parameters can be identified from the input- 
output behavior by checking the observability of the ex¬ 
tended system. Consequently, PI can be considered as a 
special observer design problem. 


2. Network Reconstruction 

When the system parameters contain information 
about the network structure, the corresponding PI prob¬ 
lem can be generalized to a network reconstruction (NR) 
problem. Consider a network whoe state variables are 
governed by a set of ODEs 

N 

^ fij J ttj(t), (^^) 

where t = 1, • • • ,7V; the coupling functions fij : M x K —>■ 
K capture the interactions between nodes: self interac¬ 
tions when i = j or pairwire interactions when i ^ j. 
The term Ui(t) € IR represents either known signals or 
control inputs that can affect node f’s state. The in¬ 
teraction matrix A = [a^] G captures the di¬ 

rected interactions between the nodes: Oji ^ 0 if node 
j directly affects node i’s dynamics. Given measured 
temporal data {xi{t),Ui{t)}fLi, Vt G [toi^i]) NR aims to 
recover some properties of the A matrix, e.g. its sign pat¬ 
tern S = [sy] = [sign(ay)] G {—1,0,1}”^", connectivity 
pattern C = [cij] = [|sij|] G {0,1}"^", adjacency pattern 
K = [kij] = [cij{l — 5ij)] G {0, !}"><" ((5j^- is the Kronecker 
delta) or in-degree sequence d = [di] = ^ 

Note that PI aims to recover the A matrix itself. 

There are three principally different NR approaches. 


which assume various levels of a priori knowledge about 
the system (Timme and Casadiego 2014). 

Driving-response. Here we try to measure and evaluate 
the collective response of a networked system to external 
perturbations or driving. As the response depends on 
both the external driving signal (which unit is perturbed, 
when and how strong is the perturbation, etc.), and the 
(unknown) structural connectivity of the network, suffi¬ 
ciently many driving-response experiments should reveal 
the entire network. This approach is relatively simple to 
implement and the required computational effort scales 
well with the system size. It has been well established for 


the reconstruction of gene regulatory networks (Gardner 


et al. 2003[ [Tegner et al.\ |2003[ |Yu[ |2010[ |Yu and Par- 


litz 


2010). Yet, this approach requires us to measure 


and drive the dynamics of all units in the system, which 
is often infeasible. The collective dynamics suitable for 
the driving-response experiments also needs to be simple 
(i.e., to exhibit a stable fixed point or periodic orbits, 
or to allow the system to be steered into such a state). 
For systems exhibiting more complex features, e.g. chaos, 
bifurcations, multi-stability, this approach is not applica¬ 
ble. If the system exhibits the same fixed point for differ¬ 
ent constant inputs (as some biological systems that have 
“perfect adaptation”), it is impossible to reconstruct the 


network using driving-response experiments (Prabakaran 


et al.. 2014). 


Copy-synchronization: This approach sets up a copy 
of the original system and updates its interaction matrix 
continuously until the copy system synchronizes its tra¬ 
jectories with the original system (Yu et al. 2006). We 
expect the hnal interaction matrix of the copy system to 
converge to that of the original system. Unfortunately, 
sufficient conditions for the convergence of this approach 
have not been fully understood and the approach is model 
dependent. Knowing the details of the coupling functions 
fij{xi,Xj) is crucial to set up the copy system. Further¬ 
more, fij{xi,Xj) needs to be Lipschitz continuous. These 
constraints significantly narrow the applicability of this 
approach. 

Direct approach: This approach relies on the 

evaluation of temporal derivatives from time series 
data (Shandilya and Timme 2011). Exploiting smooth¬ 
ness assumptions, it finds the unknown interaction ma¬ 
trix by solving an optimization problem (e.g., ii or 1 ' 2 - 
norm minimization). The rationale is as follows. If the 
time derivatives of the state variables are evaluated, and 
if the system coupling functions are also known, then the 
only remaining unknown parameters are the edge weights 
or interaction strengths Oi^’s. Repeated evaluations of 
(78) at different sufficiently closely spaced times G M 
comprise a simple and implicit restriction on the interac¬ 
tion matrix A. This approach serves as a simple starting 
strategy of NR. Yet, it has an fundamental drawback — 
there is no reason why the true interaction matrix should 
be optimal in some a priori metric. Moreover, it may suf- 
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fer from the poor evaluation of time derivatives of noisy 
time series data. 

All three approaches suffer from one common issue: 
The necessary and sufficient conditions under which they 
succeed are unknown. An important exception is the 


Modular Response Analysis method (Kholodenko et al. 


2002 Sontag 2008), which is a special driving-response 


approach, and guarantees to recover the interaction ma¬ 
trix using steady-state data collected from sufficiently 
many perturbation experiments. One drawback of this 
method is that it assumes the system is not retroactive 
(Sontag 2002 2011). Here, retroactivity manifests as 


“load” or “impedance” effects that might be hard to an¬ 
ticipate if we have no a-priori knowledge of the system 
dynamics. 

Recently, two classes of fundamental limitations of NR 
were characterized by deriving necessary (and in some 
cases sufficient) conditions to reconstruct any desired 


property of the interaction matrix (Angulo et al.. 2015). 


The first class of fundamental limitations is due to our 
uncertainty about the coupling functions fij {xi,Xj), lead¬ 
ing to a natural trade-off: the more information we want 
to reconstruct about the interaction matrix the more cer¬ 
tain we need to be about the coupling functions. For ex¬ 
ample, it is possible to reconstruct the adjacency pattern 
K without knowing exactly the coupling functions. But, 
in order to reconstruct the interaction matrix A itself, 
it is necessary to know these functions exactly. In this 
sense, if we are uncertain about the coupling functions, 
NR is easier than PI. The second class of fundamental 
limitations originates solely from uninformative tempo¬ 
ral data, i.e. {xi{t),Ui{t)}fLj^, Vt G This leads to 

a rather counterintuitive result: regardless of how much 
information we aim to reconstruct (e.g. edge weights, 
sign pattern or connectivity pattern), the measured tem¬ 
poral data needs to be equally informative. This happens 
even if we know the coupling functions exactly. Hence, 
in the sense of informativeness of the measured data, re¬ 
constructing any property of the interaction matrix is as 
difficult as reconstructing the interaction matrix itself, 
i.e. NR is as difficult as PI. A practical solution to cir¬ 
cumvent this limitation without acquiring more temporal 
data (i.e. performing more experiments, which are some¬ 
time either infeasible to too expensive), prior knowledge 
of the interaction matrix, e.g. the bounds of the edge 


weights, is extremely useful (Angulo et al. 2015). 


V. TOWARDS DESIRED FINAL STATES OR 
TRAJECTORIES 


arbitrary linear system into an arbitrary final state using 
the minimum control energy ||u(t)|pdt. For nonlin¬ 
ear dynamics we lack a ready-to-use solution, and finding 
one can be very difficult. Yet, solving such nonlinear con¬ 
trol problems has important applications from robotics 
to ecosystem management, and from cell reprogram¬ 
ming to drug discovery. For example, in robotics engi¬ 
neers frequently encounter the so-called motion- or path¬ 
planning problem, needing to decompose a desired move¬ 
ment into discrete motions that satisfy specific movement 
constraints and possibly optimize some aspect of the tra¬ 
jectory. The parallel parking problem is a typical exam¬ 
ple, requiring us to determine the sequence of motions a 
car must follow in order to parallel park into a parking 
space. 


In many cases, we are interested in steering the sys¬ 
tem towards a desired trajectory or attractor, instead of 
a desired final state. A trajectory or an orbit of a dy¬ 
namical system is a collection of points (states) in the 
state space. For example, a periodic orbit repeats itself 
in time with period T, so that x(t) = x(t -|- nT) for any 
integer n > 1. Roughly speaking, an attractor is a closed 
subset A of a dynamical system’s state space such that 
for “many” choices of initial states the system will evolve 


towards states in A (Milnor 2006). Simple attractors 


correspond to fundamental geometric objects, like points, 
lines, surfaces, spheres, toroids, manifolds, or their sim¬ 
ple combinations. Fixed (or equilibrium) point and limit 
cycle are common simple attractors. Fixed points are 
defined for mappings Xn+i = /(a^n), where a; is a fixed 
point if X = /(x), whereas equilibrium points or equi¬ 
libria are defined for flows (ODEs) x = f(x), where x 
is an equilibrium point if f(x) =0. A limit cycle is a 
periodic orbit of the dynamic system that is isolated. An 
attractor is called strange if it has a fractal structure that 
cannot be easily described as fundamental geometric ob¬ 
jects or their simple combinations. A strange attractor 
often emerges in chaotic dynamics. 


In this section we briefly review progress made in sev¬ 
eral directions with the common goal of controlling some 
dynamical systems: (a) Control of chaos, which requires 
us to transform a chaotic motion into a periodic tra¬ 


jectory using open-loop control (Hubler et al.. 1988), 


Poincare map linearization (Ott et al. 1990) or time- 


delayed feedback (Pyragas 1992). (b) Systematic design 
of compensatory perturbations of state variables that 
take advantage of the full basin of attraction of the de- 


A signihcant body of work in control theory deals with sired final state ( 

Cornelius et al. 

2013). (c) Construction 

the design of control inputs that can move the system of the attractor network ( 

Lai 

r 

014 

Wang et al. 

2014); 


from a given initial state to a desired final state in the (d) Mapping the control problem into a combinatorial op- 


state space (Sontag 1998). For linear dynamics, Equa- timization problem on the underlying networks (Eiedler 


tion (35) provides the optimal input signal to take an et al.. 2013 Mochizuki et al. 2013). 
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A. Controlling Chaos 


A deterministic dynamical system is said to be chaotic 
if its evolution is highly sensitive to its initial conditions. 
This sensitivity means that arbitrary small measurement 
errors in the initial conditions grow exponentially with 
time, destroying the long-term predictability of the sys¬ 
tem’s future state. This phenomenon, known as the 


butterfly effect, is often considered troublesome (Lorenz 


19631. Chaotic behavior commonly emerges in natural 


and engineered systems, being encountered in chemistry, 
nonlinear optics, electronics, fluid dynamics, meteorol¬ 


ogy, and biology (Strogatz 19941. 


It has been realized that well-designed control laws can 
overcome the butterfly effect, forcing chaotic systems to 


follow some desired behavior (Hubler et al.. 1988 |Ott 
et al . 1990[ Pyragas[ |1992[ |Sass and Toroczkai[ |1996f 
Toroczkai, 19941. Next, we review several key methods 


devised for the control of chaotic systems from the control 


theoretical perspective ( 

Boccaletti et al. 

2000 

Chen and 

Dong 1998 

Fradkov and Evans 

2005 

I. 


1. Open-loop Control 


Since the late 1980s, a series of methods have emerged 
to manipulate chaotic systems towards a desired “goal 
dynamics” g(t) (Hubler et al. 1988). Consider a con¬ 
trolled system 


X = F(x) -t Bu{t) (79) 

where x C is the state vector, u(t) G is the con¬ 
trol input. In contrast with the network-based problems 
discussed earlier, here we assume that all state variables 
are controlled (M = N) and detB 0. The goal is 
to design u{t) so that x(f) converges to a desired tra¬ 
jectory g(t), i.e., |x(t) — g(t)| —>• 0 as t —>■ oo. We can 
use open-loop control for this purpose, using the control 
input called the Hubler action, 


Note that the method (79l-(80) is not tailored to 
chaotic systems, but potentially works for any nonlin¬ 
ear system. It has several disadvantages, though: (i) 
the open-loop control (801 requires a priori knowledge 
of the dynamics, which is often not precisely known for 
complex systems; (ii) the applied controls are not always 
small, requiring high control energy; (iii) the convergence 
of |x(t) — g(t)| —>• 0 for t —oo depends on the detailed 
functional form of F(x) and the initial condition x(0), 
hence this method is not guranteed to work for arbitrary 
systems. 


2. Linearization of the Poincare map: OGY method 


The OGY method proposed by Ott, Grebogi and Yorke 


(Ott et al. 1990) exploits the observation that typically 
an infinite number of unstable periodic orbits (UPOs) 
are embedded in a chaotic attractor (Fig. 29). Therefore 
we can obtain a desired periodic motion by making only 
small perturbations to an accessible system parameter. 

The OGY method can be summarized as follows: First, 
we determine and examine some of the low-period UPOs 
embedded in the chaotic attractor. Second, we choose a 
desired UPO. Finally, we design small time-dependent 
parameter perturbations to stabilize this pre-existing 
UPO. 


This method is not only very general and practical, but 
also suggests that in some systems the presence of chaotic 
behavior can be an advantage for control. Indeed, if the 
attractor of a system is not chaotic but has a stable pe¬ 
riodic orbit (SPO), then small parameter perturbations 
can only slightly change the existing orbit. Therefore, 
given that any one of the infinite number of UPOs can be 
stabilized, we can always choose the UPO that achieves 
the best system performance. Hence, chaotic behavior 
offers us a diverse and rich landscape for the desired dy¬ 
namic behavior of the system. 

To demonstrate this method, let us consider a nonlin¬ 
ear continuous-time dynamical system 


u(t) = B-Mg(t)-F(g(f))], (80) 


X = f(x, u) 


(82) 


which ensures that x(t) = g(t) is a solution of the con¬ 
trolled system. In this case, the error e(t) = x(t) — g(t) 
satisfies 


e{t) =F{e{t) + g{t))-F{g{t)), (81) 

which can be linearized as e{t) = A{t)e{t), where A(t) = 
|x=g(t)- If the linearized system is uniformly asymp¬ 
totically stable, i.e., its equilibrium point e* = 0 is stable 
for all t > 0, then the error e{t) converges to zero, and 
x(t) converges to the desired trajectory g(t). We call 
the regions of the state space from which the controlled 
orbits converge to the goal trajectory g{t) entrainment 
regions. 


where x G R^ is the state vector and u S M repre¬ 
sents a tunable parameter, which can be considered as 
a control input. Our task is to reach a desired trajec¬ 
tory x*(t) that satisfies (82) with it = 0. To achieve that 
we first construct a surface S, called a Poincare section, 
which passes through the point Xq = x* (0) transversally 
to the trajectory x*(t) (see Fig. 30). Consider a map 
X i~> F(x, It), where F(x, it) is the point of first return to 
the Poincare section of the solution of (821 that begins at 
the point x and was obtained for the constant input it. 
Since we can integrate (82) forward in time from x, the 
map X F(x,it), called the Poincare map, must exist. 
Note that even though we may not be able to write down 
the map F explicitly, the knowledge that it exists is still 
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FIG. 29 Chaotic behavior in a nonlinear electronic circnit. 
The vertical axis measures the voltage drop V (t) across a 5017 
resistor, being proportional to the current in the circuit. The 
system ergodically visits the unstable periodic orbits (UPOs) 
embedded in t he chaotic attractor . The plot shows three such 
UPOs. After (Sukow et al. 19971. 



FIG. 30 (Color online) Poincare map. In a continuous dy¬ 
namical system the Poincare map is the intersection of a peri¬ 
odic orbit in the state space with a certain lower-dimensional 
subspace, called the Poincare section S, transversal to the 
flow of the system. In the Poincare section S, the Poincare 
map X I— >■ F(x,u) projects point x onto point F(x,u), i.e., 
Xk = F(xfc_i,Ufc_i), Xfc+i = F(xfc,UA,), • ■ •. 


where zj, = x^ — Xq, A = ^ is the Jacobian matrix, 

OX \ Xq 

and B = I is a column vector. 

OU Ixo 


To stabilize the linear system (84) and hence steer the 


original system to a desired periodic orbit that passes 
through Xq, the OGY method employs a linear state feed¬ 
back control law 


Uk = 


Czk if |zfe| < (5 
0 otherwise 


(85) 


where J > 0 is a sufficiently small parameter. Note that 
the control is only applied in some neighborhood of the 
desired trajectory, which ensures the smallness of the con¬ 
trol action. This piecewise-constant small action control 
is a key feature of the OGY method. To guarantee the 
efficiency of the method, the matrix C must be chosen so 
that in the linear closed-loop system Zk+i = (A+BC) z^, 
the norm |(A + BC)z| < p|z| decreases, where p < 1. 

Extensive numerical simulations have corroborated the 
practical utility of the OGY method. Furthermore, the 
OGY method was proven to be effective in experimen¬ 
tal systems as well, allowing the stabilization of unstable 
periodic orbits in a chaotically oscillating magnetoelastic 
ribbon, a driven diode circuit, a multimode laser with an 
intracavity crystal, a thermal convection loop, and the 


Belousov-Zhabotinsky reaction (Boccaletti et al. 20001. 


Slow convergence was often reported, a price we must 
pay to achieve global stabilization of a nonlinear system 


with small control action (Fradkov and Evans 20051. 


The advantage of the OGY method is that it does not 
require prior knowledge of the system’s dynamics. In¬ 
stead, we just rely on the system’s behavior to learn the 
necessary small perturbation to nudge it towards a de¬ 
sired trajectory. This is similar to the balancing of a stick 
on our palm, which can be achieved without knowing 
Newton’s second law of motion and the stick’s detailed 
equation of motion. Indeed, in the OGY method, both A 


and B in (841 can be extracted purely from observations 


of the trajectory on the chaotic attractor (Shinbrot et al. 


1993). Finally, the OGY method can be extended to 


arbitrarily high dimensional systems, without assuming 


knowledge of the underlying dynamics (Auerbach et al. 


1992). 


useful (Shinbrot et al. 1993). By considering a sequence 


of such maps, we get a discrete system 

Xfe+i = F(xfc,Ufe), 


(83) 


where x^ = x(tfc), tk is the time of the fc-th intersection 
of the Poincare section S, and Uk is the value of control 
u{t) over the interval between tk and tk+i. 

A key step in the OGY method is to linearize the dis¬ 


crete system (83) as 


3. Time-delayed feedback: Pyragas method 

The Pyragas method employs continuous feedback to 
synchronize the current state of a system with a time- 
delayed version of itself, offering an alternative approach 
to stabilizing a desired UPO embedded in a chaotic at¬ 
tractor ( {Pyragas' 1992). Consider the nonlinear system 
(82). If it has a desired UPO T = {x*(t)} with period T 
for u = 0, then we can use the feedback control 


(84) 


u{t) = K [x(t) — x(t — r) 


Z/c + l AZ/;; + BUfc, 


( 86 ) 
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FIG. 31 (Color online) Controlling chaos. The image shows 
the use of the Ott-Grebogi-Yorke (OGY) method to control 
the chaotic behavior in the Henon map Xn+i = p + 0.3Yn — 
Y2 ,Y„+i = X„ where the parameter p is set to po = 1.4. 
(a) The Henon attractor contains period-1 point A*, which 
is revisited in each map iteration, period-2 points Bi and 
B 2 , which are revisited every other map iteration, i.e., Bi —>■ 
B 2 —>■ Hi —>■ 132 , and period-4 points Ci,C2,C3 and 

C4, which are cycled through every four map iterations, (b) 
The result of stabilizing the periodic orbit A* of the Henon 
attractor by tuning p by less than 1% around po- The arrow 
indicates the time step at which the small perturbation is 
initiated. For the first 86 iterations, the trajectory moves 
chaotically on the attractor, never falling within the desired 
small region about A*. On the 87th iteration, following the 
application of the control perturbation, the s tate falls within 
the desired region, and is held near A*. After | Shinbrot et al. 


19931. 


where K is the feedback gain and r is the delay time, to 
stabilize the desired UPO. If r = T and the solution x(f) 


of the closed-loop system (82 86) begins on the UPO, 


then it remains on the UPO for all t > 0. Surprisingly, 
x(t) can converge to the UPO even if initially is not on 
the UPO, i.e., x(0) ^ P. 

Considering that not all the state variables are exper¬ 


imentally accessible, we can rewrite (86) as 


L{t) = K [y{t) - y{t - T)] 


(87) 


for a desired UPO of period T. Here y(t) = h(a:(t)) € M 
is an experimentally accessible output signal. The advan¬ 


tage of the time-delayed feedback control law (87) is that 


it does not require rapid switching or sampling, nor does 
it require a reference signal corresponding to the desired 
UPO. Unfortunately, the domain of system parameters 


over which control can be achieved via (87) is limited. 


Furthermore, the method fails for highly unstable orbits. 
Note, however, that an extended variant of the Pyragas 
method, using a control law whose form is closely related 
to the amplitude of light reflected from a Fabry-Perot in¬ 


terferometer can stabilize highly unstable orbits (Socolar 


et al. 1994). 


Despite the simple form ( 86|87 ) of the control signal, 
the analytical study of the closed-loop system is chal¬ 
lenging. Indeed, while there are extensive numerical and 


experimental results pertaining to the properties and ap¬ 
plication of the Pyragas method, the sufficient conditions 


that guarantee its applicability remain unknown (Frad¬ 


kov and Evans 2005). 


Note that similar to the Pyragas method, a geometric 


method of stabilizing UPOs (Sass and Toroczkai, 1996 


Toroczkai, 1994) also uses time-delays. This method is 


based on some observations about the geometry of the 
linearized dynamics around these orbits in the phase 
space. It does not require explicit knowledge of the 
dynamics (which is similar to the OGY method), but 
only experimentally accessible state information within 
a short period of the system’s immediate past. More 
specifially, it requires a rough location of the UPO and a 
single parameter easily computed from four data points. 
This geometric method does not have the problems of the 
Pyragas method in stabilizing UPOs. The drawback of 
this geometric method is that it has only been formulated 
for 2D maps and 3D flows. 


B. Compensatory Perturbations of State Variables 


The control tools described above were mainly de¬ 
signed for low dimensional dynamical systems with a 
simple structure. Most complex systems are high¬ 
dimensional, however, consisting of a network of com¬ 
ponents connected by nonlinear interactions. We need, 
therefore, tools to bring a networked system to a desired 
target state. A recently proposed method can work even 
when the target state is not directly accessible due to 
certain constraints (Cornelius et al. 2013). The basic in¬ 


sight of the approach is that each desirable state has a 
“basin of attraction”, representing a region of initial con¬ 
ditions whose trajectories converge to it. For a system 
that is in an undesired state, we need to identify pertur¬ 
bations to the state variables that can bring the system to 
the basin of attraction of the desired target state. Once 
there, the system will evolve spontaneously to the target 
state. Assume that a physically admissible perturbation 
fulfills some constraints that can be represented by vector 
expressions of the form 


g(xo,Xo) < 0 and h(xo,Xo) = 0, (88) 

where the equality and inequality apply to each compo¬ 
nent individually. To iteratively identify compensatory 
perturbations we use the following procedure: Given the 
current initial state of the network, Xq, we integrate the 
system’s dynamics over a time window to < t < to + T 
to identify the time when the orbit is closest to the 
target, U = argminjx* — x(t)|. We then integrate 
the variational equation up to U to obtain the corre¬ 
sponding variational matrix, M(tc), which maps a small 
change i5xo in the initial state of the network to a change 
(5x(tc) in the resulting perturbed orbit at U according 
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FIG. 32 (Color online) Steering a network using compen¬ 
satory perturbations of state variables, (a) The control set 
(shown in yellow) is a set of nodes that are accessible to com¬ 
pensatory perturbations, (b) In the absence of control, the 
network is in an initial state xq and evolves to an undesir¬ 
able equilibrium Xu (red curve). By perturbing the initial 
state (orange arrow), the network reaches a new state Xg, 
which evolves to the desired target state (blue curve), (c) 
Typically, the compensatory perturbations must obey some 
constraints. In this example, we can only perturb three out 
of N dimensions, corresponding to the three-node control set 
(shown in yellow), and the state variable of each control node 
can only be reduced. These constraints form a cube (grey 
volume) within the three-dimensional subspace of the control 
nodes. The network can be steered to the target state if and 
only if the corresponding slice of the target’s basin of attrac¬ 
tion (blue volume) intersects this cube, (d) Along each orbit 
there is a point that is closest to the target state. We seek to 
identify a perturbation (magenta arrow) to the initial condi¬ 
tion that brings the closest point closer to the target (green 
arrow), (e) This process is repeated (dashed curves), until we 
identify a perturbed state x'q that is in the attraction basin 
of the target state, hence the system automatically evolves to 
the target state. This results in a compensatory pertur bation 
xo —> Xq (orange arrow). After (Cornelius et al. 20131. 


to 6 x{tc) = M(tc) • ^xq. This mapping is used to se¬ 
lect an incremental perturbation Sxq that minimizes the 
distance between the perturbed orbit and the target at 
time tc, subject to the constraints (881 and additional 
constraints on ^Xq to ensure the validity of the varia¬ 
tional approximation. 

The selection of Jxq is performed via a nonlinear op¬ 
timization that can be efficiently solved using sequential 
quadratic programming. The initial condition is then up¬ 
dated Xq —> Xg -|-<5xo, and we test whether the new initial 
state lies in the target’s basin of attraction by integrating 
the system dynamics over a long time r. If the system’s 
orbit reaches a small ball of radius k around x* within 
this time, we declare success and recognize Xg — Xg as a 
compensatory perturbation (for the updated Xg). If not, 
we calculate the time of closest approach of the new orbit 
and repeat the procedure. 

Similar to the open-loop control of chaos discussed in 
Sec. V.A.l the approach based on compensatory pertur¬ 


bation potentially works for any nonlinear system. It has 
been successfully applied to the mitigation of cascading 
failures in a power grid and the identification of drug 


targets in a cancer signaling network (Cornelius et al. 


2013). Yet, the approach requires a priori knowledge of 


the detailed model describing the dynamics of the system 
we wish to control, a piece of knowledge we often lack in 
complex systems. With an imperfect model, a compen¬ 
satory perturbation may steer the system into a different 
basin of attraction than the desired one. Studying the 
dependence of the success rate of this approach on the 
parameter uncertainty and system noise remains an ana¬ 
lytically challenging issue. Moreover, it is unclear how to 
choose the optimal control set consisting of one or more 
nodes accessible to compensatory perturbations so that 
some control objectives, like the number of control nodes 
or the amount of control energy, are minimized. 


C. Small Perturbations to System Parameters 


The control tool described above perturbs the state 
variables of a networked system. In analogy with the 


OGY method (Ott et al.. 1990), we can also control a 


networked system via small perturbations to its param¬ 
eters. Note that networked systems are typically high¬ 
dimensional, to which the control methodologies devel¬ 
oped for chaotic system do not immediately apply. Yet, 
we can control complex networked systems via perturba¬ 


tions to the system parameters (Lai 2014), an approach 


complementary to the approaches based on perturbations 
of the state variables. The key step of this approach 
is to choose a set of experimentally adjustable parame¬ 
ters and determine whether small perturbations to these 
parameters can steer the system towards the desired at¬ 


tractor (Lai, 1996 2014). Depending on the physical con¬ 


straints the control parameters obey, the directed control 
path from the undesired attractor to the desired attractor 
can either be via a direct connection or via intermediate 
attractors along the control path. If there are no feasible 
control paths reaching the desired attractor, then we can 
not steer the system to that attractor, hence control is 
not possible. 

Considering each attractor as a node, and the con¬ 
trol paths as directed edges between them, we can con¬ 
struct an “attractor network”, whose properties deter¬ 
mine the controllability of the original dynamic net¬ 


work (Lai 2014). For a given nonlinear system, the at¬ 
tractor network can be constructed as follows. First, we 
identify all possible attractors and choose a set of sys¬ 
tem parameters that can be experimentally perturbed. 
Second, we set the system into a specific attractor a, and 
determine the set of attractors into which the system can 
evolve from the original attractor a with a reasonable 
combination of the adjustable parameters. This effec¬ 
tively draws a link from attractor a to all other attractors 
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FIG. 33 (Color online) Epigenetic state network (ESN), (a) 
On the epigenetic landscape, a minimal energy path connects 
two neighboring attractors through an unstable transition 
point (first-order saddle point). The landscape can be rep¬ 
resented by a network, where nodes are attractors or basins 
of attraction and edges are minimal energy paths connecting 
the neighboring attractors, (b) The vector held of a mutu¬ 
ally inhibitive two-gene circuit (inset). Nodes SI, S2 and S3 
are hxed-point attractors. The pie diagram of each attrac¬ 
tor represents the expression pattern of the two genes. The 
hrst-order saddle points (green diamond) are surrounded by 
forward and backward optimal paths (dark blue) connecting 
two neighboring attractors, (c) The ESN constructed from 


(a) by connecting neighboring attractors. After (Wang et al. 


20141. 


reachable by feasible parameter perturbations. Finally, 
we repeat this procedure for all attractors, obtaining the 


attractor network (Lai 20141. 


To illustrate the construction of such an attractor net¬ 
work, consider the epigenetic state network (ESN) that 
describes the phenotypic transitions on the epigenetic 


landscape of a cell (Fig. 331. In the epigenetic landscape, 
two neighboring fixed-point attractors, corresponding to 
stable cell phenotypes, are connected by a minimal en¬ 
ergy path through an unstable transition point (first- 


order saddle point) (Wang et oL, 2011 20141. The num¬ 


ber of fixed points (nodes) and saddle points (edges) 
grows exponentially with the number of genes (dimen¬ 
sionality). We can rely on a conditional root-finding 


algorithm (Wang et al.\ 20141 to construct this epige¬ 


netic state network (ESN). The obtained ESN captures 
the global architecture of stable cell phenotypes, help¬ 
ing us translate the metaphorical Waddington epigenetic 


landscape concept (Slack, 2002 Waddington and Kacser 


19571 into a mathematical framework of cell phenotypic 


transitions. 


D. Dynamics and Control at Feedback Vertex Sets 


For regulatory networks described as a digraph of de¬ 
pendencies, it has been recently shown that open-loop 
control applied to a feedback vertex set (FVS) will force 
the remaining network to stably follow the desired tra¬ 


jectories (Fiedler et aL, 2013 Mochizuki et al.\ 2013). An 


FVS is a subset of nodes in the absence of which the di¬ 
graph becomes acyclic, i.e., contains no directed cycles 


(Fig. 34). Unlike the approaches discussed in Sec. V.B 


and Sec. |V.C[ this approach has rigorous analytical un¬ 
derpinnings. 

Consider a general non-autonomous nonlinear net¬ 
worked system 


Xi = Fi{t,Xi,xi.) 


(89) 


where i = 1, • • • , and li denotes the set of upstream 
or incoming neighbors of node 1, i.e., j S Ti if there is 
directed edge (j i) in the network. The correspond¬ 
ing network is often called the system digraph, which is a 
transpose of the inference diagram introduced in observ¬ 


ability (see Sec. IV.B). 


An open-loop control applied to the nodes of an FVS 
will completely control the dynamics of those nodes and 
hence effectively remove all the incoming links to them. 
Consequently, those nodes will not be influenced by other 
nodes. They will, however, continue to influence other 
nodes and drive the whole system to a desired attractor. 
Consider, for example, the gene regulatory network of cir¬ 


cadian rhythms in mice, consisting of 21 nodes (Fig. 35 1 ). 
In general there can be multiple minimal FVS’s for a 
given digraph. One such minimal FVS of size seven, 
i.e. F = {PERI, PER2, CRYl, CRY2, RORc, CLK, 


BMALl}, is highlighted in red in Pig. 35 1 . The associ¬ 


ated dynamical system can be described by a set of ODEs 
involving 21 variables and hundreds of parameters. Un¬ 
der a particular choice of parameters, the system has 
several invariant sets: (i) two stable periodic oscillations 
(PI and P2); (ii) one unstable periodic oscillation (UP); 
and (hi) one unstable stationary point (USS) (Fig.[3^,c). 
Let us aim to steer the system from PI to P2. To achieve 
this, we first need to calculate the time tracks of the seven 
FVS nodes on the desired invariant set P2, denoted as 
xf^,i G F, which can be done by numerically integrat¬ 
ing the ODEs. Then we prescribe the time tracks of the 
seven nodes in F to follow their desired values xf^. This 
way, we effectively remove any influence from the other 
14 nodes to the nodes in F. The dynamics of the remain¬ 
ing 14 nodes Xi,i ^ F, are determined by the remaining 
14 ODEs of the system, where the initial state of these 
remaining nodes is chosen to coincide with an arbitrary 
point on the PI trajectory. As shown in Fig. [35jl, the 
trajectories of the remaining 14 nodes deviate from the 
original stable periodic orbit PI and quickly converge to 
the competing orbit P2. The whole system eventually 
displays periodic oscillation on the P2 orbit. 
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In the above example, the identified FVS is a minimal 
one, i.e., any subset of T is not an FVS. Yet, a mini¬ 
mal FVS is not guaranteed to be the minimum one that 
contains the least number of nodes. Naturally, it will 
be more desirable to identify and control the nodes in 
the minimum FVS. Unfortunately, finding the minimum 


FVS of a general digraph is an NP-hard problem (Karp 


19721. 


This FVS-based open-loop control can be applied to a 
wide range of nonlinear dynamical systems. It requires 
only a few conditions (e.g. continuous, dissipative and 
decaying) on the nonlinear functions F) that are very 


mild and satisfied by many real systems (Fiedler et al. 


20131. 


For systems associated with a digraph G{V, E), we can 
rigorously prove that clumping the dynamics of a sub¬ 
set of nodes S' C U will control the rest of the network 
towards the desired attractor for all choices of nonlin¬ 
earities Fi that satisfy the above-mentioned conditions if 
and only if S is an FVS in G (Fiedler et al. 20131. Yet, 
there do exist specific systems (with certain nonlinearity 
Fi) where clumping a reduced FVS (i.e. removing one 
or more nodes from an FVS) is sufficient to control the 
system to a desired attractor. In other words, for a spe¬ 
cific system, clumping an FVS might be not necessary. 
It would be a natural starting point, though. 

Note that to apply the two approaches discussed in 
the previous subsections, namely the compensatory per¬ 
turbations of state variables (Sec. V.Bl, and attractor 
network based on small perturbations of system parame¬ 
ters (Sec. V.C I, we need a detailed knowledge of the sys¬ 
tem dynamics, including all system parameters. In many 
cases, we lack such a piece of knowledge. In contrast, to 
apply the FVS-based open-loop control (Sec. V.D), we 
just need the trajectories of FVS nodes on the desired 
attractors. We do not have to know full dynamics, nor 
the exact parameter values. We just need to assure a few 
mild conditions on the nonlinear functions Fi are satis¬ 
fied. 


VI. CONTROLLING COLLECTIVE BEHAVIOR 


Dynamical agents interacting through complex net¬ 
works can display a wide range of collective behavior, 
from synchronization to flocking among many interacting 
agents. In particular the study of network-mediated syn¬ 
chronization has a long history, with applications from 
biology to neuroscience, engineering, computer science. 


economy and social sciences (Arenas et al. , 20081. Flock¬ 


ing has also gained significant attention in the past two 
decades, capturing phenomena from the coordinated mo¬ 
tion of birds or fish to self-organized networks of mo¬ 
bile agents. Applications range from massive distributed 
sensing using mobile sensor networks to the self-assembly 
of connected mobile networks, and military missions such 
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FIG. 34 Feedback vertex set (FVS). This figure show exam¬ 
ples of FVSs in directed graphs, whose removal render the 
graphs acyclic. The gray vertices represent a choice of a min¬ 
imal FVS in each panel (a)-(e). Controlling the dynamics 
of the nodes in an FVS allows us to switch the dynamics of 
the whole system from one attractor to some other attractor. 
After | |Mochizuki et a^|2013| ). 


as reconnaissance, surveillance, and combat using co¬ 


operative unmanned aerial vehicles (Olfati-Saber et al. 


2007). These problems pose, however, a number of 


fundamental questions pertaining to the control of self- 
organized networks. 


If we aim to achieve a desired collective behavior, it 
is often infeasible to directly control all nodes of a large 
network. This difficulty is partially alleviated by the no¬ 
tion of pinning control (Wang and Chen 2002a|b ), which 
relies heavily on feedback processes. In pinning control 
a feedback control input is applied to a small subset of 
nodes called pinned nodes, which propagates to the rest 
of the network through the edges. The design and im¬ 
plementation of feedback control must take into account 
both the individual dynamics of the components and the 
network topology. Conceptually, pinning control is sim¬ 
ilar to the minimum controllability problem of a linear 
system discussed in Sec.|^ The key difference is that, in¬ 
stead of fully controlling a system, pinning control aims 
to control only the system’s collective behavior, like syn¬ 
chronization or flocking. Pinning control has been ex¬ 
tensively applied to the synchronization of coupled oscil¬ 


lators and flocking of interacting agents (Bernardo and 


DeLellisl 2014| Chen and Duan 2008 Li et al.\ 2004 For- 

firi and di Bernardol 2008 Sorrentino et al.\ 2007a 

Wang 

and Chen 2002a|bl Yu et al. 2009a 20131 Zou anc 

Chenl 


2008). 


In this section we review some fundamental results on 
controlling the collective behavior of complex networked 
systems. We pay particular attention to the pinning con¬ 
trol of synchronization and flocking. Synchronization 
of coupled oscillators is typically studied on fixed net¬ 
work topology. We build on the master stability formal¬ 
ism to explore pinning synchronization, focusing on local 
and global stability conditions and adaptive strategies. 
Flocking of multi-agent systems are typically associated 
with switching or time-varying network topology, because 
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FIG. 35 (Color online) Controlling a system through its feedback vertex set (FVS). (a) A regulatory network with 21 variables 
describes the mammalian circadian rhythms in mice ( |Mirsky et al. 20091. A minimal FVS of seven elements, denoted as T, is 
highlighted by red circles, (b) Trajectories of two stable periodic orbits, periodl (PI, dotted and broken curve) and period2 (P2, 
dotted curve), one unstable periodic orbit (UP, broken curve) and one unstable stationary state (USS, solid line), represented 
by time tracks of the variable Per2. (c) Trajectories of the same solutions in the phase plane of the two variables Perl and 
Per2, which are not in the FVS. (d-g) Numerical trajectories of successful open loop controls of circadian rhythms via the full 
feedback vertex set T. Zooms into P2, UP, and USS are shown as top-right insets. The resulting trajectory of the control 
experiment is always the red solid curve, (d) “From PI to P2”. The stable cycles PI and P2 are shown by gray solid curves, 
(e) “From P2 to PI”. Gray soli d: PI and P2. (f ) “From PI to UP”. Gray solid: PI and UP. (g) “From PI to USS”. Gray 
solid: PI, open dot: USS. After ( [Mochizuki et a/!||2013[ ). 


the agents, like robots, vehicles or animals, are often mo¬ 
bile. To illustrate this we discuss the Vicsek model of 
flocking behavior, emphasizing its control theoretical in¬ 
terpretation. Finally, we review key protocols that can 
induce flocking in multi-agent systems. 


A. Synchronization of coupled oscillators 

Consider a static network of N identical nodes (oscil¬ 
lators) with nearest-neighbor coupling: 

N 

Xi = f(xi) + a'^a^jWij[h{xj) - h(xi)] 

N 

= i{xi)-a'^g^jh{xj), (90) 

where x, € is the d-dimensional state vector of the 
fth node, f(xi) : —)■ determines the individual 

dynamics of each node, a is the coupling strength, also 


called the coupling gain, A = (arj) is the N x N ad¬ 
jacency matrix of the network, Wij > 0 is the weight 
of link {i,j). The output function h(x) : i is 

used to couple the oscillators and is identical for all os¬ 
cillators. For example, if we use h(x) = (a:, 0,0)"’’ for a 
three-dimensional oscillator, like the Lorenz or Rossler 
oscillator, it means that the oscillators are coupled only 
through their x-components. In general, h(x) can be 
any linear or nonlinear mapping of the state vector x. 
G = {gij) is the N X N coupling matrix of the network 

i9ij = -aijWij for i 7 ^ j and gu = - 5 u)- If 

Wij = 1 for all links, G is the Laplacian matrix L of the 
network. Note that G is not necessarily symmetric. 


The system (90) is synchronized when the trajectories 


of all nodes converge to a common trajectory, i.e. 
lim |jx,(t) -Xj(t)|| = 0 

t —>-oo 


(91) 


for all i,j = I,-- - ,N. Such synchronization behavior 
describes a continuous system that has a uniform move¬ 
ment, used to model synchronized neurons, lasers and 
electronic circuits (iPecora and Carroll 1998). 
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Due to the diffusive coupling, the completely synchro¬ 
nized state xi(t) = X 2 (t) = ••• = K]\[{t) = s{t) is a 


natural solution of Eq. (90). This also defines a linear 


invariant manifold, called the synchronization manifold, 
where all the oscillators evolve synchronously as s = f (s). 
Note that s{t) may be an equilibrium point, a periodic 
orbit, or even a chaotic solution. 

Despite the fact that the completely synchronized state 


is a natural solution of Eq. (90), it may not emerge spon¬ 


taneously. For example, if the coupling gain a is close 
to zero, the oscillators tend to behave independently. If 
the coupling gain cr is too strong, the oscillators may not 
synchronize either. Our goal is to identify the conditions 


under which the system (901 can synchronize. A broad 


spectrum of methods allows us to address this question 


(|Barahona and Pecora 2002 Belykh et al. 

2004b Chen 

2007 

IPecora and Carroll 

1998| Russo anc 

Di Bernardo 

2009 

Wu and Chua 

1994 

). The best-known method, dis- 


cussed next, is based on the calculation of the eigenvalues 
of the coupling matrix. 


1. Master stability formalism and beyond 


Consider the stability of the synchronization manifold 
in the presence of a small perturbation Xi(t) = s(t) -|- 
Sxi{t). By expanding f(xi) and h(xi) to the first order 
of 5xi, we obtain a linear variational equation for 6xi(t), 

N 

Sxi = J{s)5x^ - a'^gij£{s)Sxj, (92) 

i=i 

with Jacobian matrices ff{s) = |x=s and £(s) = 
|x=s- Let dX = [Jxi, • • • , dx^r]"'". Then formally we 

have 


JX = [I (g) J{s) - ctG O £:(s)] jx 


(93) 


where I is the N x N identity matrix and 0 is the Kro- 
necker product (a.k.a. matrix direct product). 

The key idea of the master stability formalism is that 
we need to consider only variations that are transverse 
to the synchronization manifold, as variations along s[t) 
leave the system in the synchronized state (Barahona 


and Pecora 2002 Pecora and Carroll, 1998). If these 


transverse variations damp out, then the synchronization 
manifold is stable. To separate out the transverse varia¬ 
tions, we can project JX into the eigenspace spanned by 
the eigenvectors of the coupling matrix G, i.e., JX = 
(P 0 Id) H with G P = G = Diag(Ai, A 2 , • • • ,Xn)- 
Then we have 


H= I 0 J'(s) — (tG 0 £(s) H, 


(94) 


which results in a block diagonalized variational equa¬ 
tion with N blocks, corresponding to N decoupled eigen- 
modes. Each block has the form 


L = [^(s) - aA£:(s)] (95) 


where is the eigenmode associated with the eigenvalue 


Xi of G. Note that in deriving (94) we have implicitly 


assnmed that the coupling matrix G is diagonalizable, 
which is always true for symmetric G. Thus each eigen¬ 
mode of the perturbation is decoupled from the others, 
and will damp out independently and simultaneously. If 
G is not diagonalizable, we can transform G into the Jor¬ 
dan canonical form. In this case, some eigenmodes of the 


perturbation may suffer from a long transient (Nishikawa 
and Motterf 2006). 


We can order the eigenvalues of G such that 0 = Ai < 
ReA 2 < • • ■ < ReA^f. Because the row sum of G is zero, 
the minimal eigenvalue Ai is always zero with the corre¬ 
sponding eigenvector ei = (1,1,..., 1)"^. Hence the first 
eigenmode = > 7 ( 8)11 corresponds to the perturba¬ 
tion parallel to the synchronization manifold. Due to the 
Gerschgorin Circle Theorem, all other eigenvalues must 
have non-negative real parts. The corresponding (N — 1) 
eigenmodes are transverse to the synchronization man¬ 
ifold and must decay to have a stable synchronization 
manifold. 


The form of each block in (951 is the same up to the 


scalar multiplier aXi. This leads to the variational equa¬ 
tion, called the master stability equation, 


i = [J + 


(96) 


For small^wehave ||£(t)|| ~ exp[A(a,/3)t], which decays 
exponentially if the maximum Lyapunov characteristic 
exponent A(a,/3) < 0. Consequently, A{a,/3) is called 
the master stability funetion (MSF). Given a coupling 
strength a, the sign of the MSF in the point aXi in the 
complex plane reveals the stability of that eigenmode. If 
all eigenmodes are stable (i.e. A(CTAi) < 0 for all Fs), then 
the synchronization manifold is stable at that coupling 
strength. Note that since the master stability formal¬ 
ism only assesses the linear stability of the synchronized 
state, it only yields the necessary, but not the sufficient 
condition for synchronization. 

For undirected and unweighted networks, the coupling 
matrix G is symmetric and all its eigenvalues are real, 
simplifying the stability analysis. In this case, depending 
on >7 and £, the MSF A(a) can be classified as follows: 

(i) Bounded: A(q;) < 0 for ai < a < a 2 . This usu¬ 
ally happens when h(x) 7 x;. The linear stability of the 
synchronized manifold requires that ai < 17 X 2 < • • ■ < 
ctAjv < ci 2 - This condition can be only fulfilled for a 
when the eigenratio R satisfies 


_ Aat Q!2 

K = < —. 

A2 Q!i 


(97) 


The beauty of this inequality comes from the fact that 
its r.h.s. depends only on the dynamics while its l.h.s. 
depends only on the network structure. If R > a 2 loi\, 
for any a the synchronization manifold is unstable, indi¬ 
cating that it is impossible to synchronize the network. 
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If i? < ailci\, the synchronization manifold is stable for 
CTmin = OixjXi < (T < tJmax = 0 : 2 /Xn- The synchroTiiz- 
ability of the network can be quantified by the relative 
interval crmax/cmin = a 2 /{aiR). A network is more syn- 
chronizable for higher (Tmax/cmin (or smaller R). 

(ii) Unbounded: A(q:) < 0 for q: > ai. The stability 
criteria of the synchronized manifold is ai < aX 2 < ■ • • < 
ctAjv, which is true if 

CT > CTmin = ai/A2. (98) 


The larger is A 2 the smaller is the synchronization thresh¬ 
old CTmin, hence the more synchronizable is the network. 


Inequalities (97) and (981 demonstrate that the MSF 


framework provides an objective criteria (i? or A 2 ) to as¬ 
sess the synchronizability of complex networks based on 
the spectrum of the coupling matrix G only, without re¬ 
ferring to specific oscillators and output functions. The 
MSF framework allows us to address the impact of the 
network topology and edge weights on synchronizabil¬ 


ity (Arenas et al. 2008). Consequently, there have been 


numerous numerical attempts to relate the spectral prop¬ 
erties of network models to a single structural character¬ 
istic of networks, like mean degree, degree heterogeneity, 
path lengths, clustering coefficient, degree-degree corre¬ 


lations, etc. (Arenas et al. 20081. The outcome of these 


analyses is occasionally confusing, because in a networked 
environment it is usually impossible to isolate a single 
structural characteristic while keeping the others fixed. 
Overall, several network characteristics can influence syn¬ 
chronizability, but none of them is an exclusive factor in 
the observed dependencies. 

The fundamental limitation of MSF is that it only as¬ 
sesses the linear or loeal stability of the synchronized 
state, which is a necessary^ but not a suffieient condition 


for synchronization (Arenas et al. 2008). To obtain a 


sufficient condition, one can use global stability analysis. 


like Lyapunov’s direct method (jBelykh et al ., 2006 2005 

2004a|b| Chml 2006112007 2008 Li et al. 2009 Wu and 

Chua 1994 

1995a|b|c) or contraction theory (Aminzare 

and Sontag 

2015| Li et al.\ 2007 Lohmilier and Slotine 

1998 

Pham 

and Slotinel 20071 Russo and Di Bernardo 

2009 

Tabareau et al.. 2010 Wang and Slotine 2005). 


2. Pinning synclironizability 

If a network of coupled oscillators can not synchro¬ 
nize spontaneously, we can design controllers that, ap¬ 
plied to a subset of pinned nodes C, help synchronize 
the network. Hence the pinned nodes behave like lead- 


ers ( 

Li et al. 

2004 

Wang and Slotine 

2005 

2006 

Wang 

and Chen 

2002a 

), forcing the remaining follower nodes 


to synchronize. This procedure, known as pinning syn- 
ehronization, is fundamentally different from spontaneous 
synchronization of coupled oscillators, where we don’t 


specify the synchronized trajectory s(t), hence the sys¬ 
tem “self-organizes” into the synchronized trajectory un¬ 
der appropriate conditions. In pinning synchronization, 
we choose the desired trajectory s(t), aiming to achieve 
some desired control objective, and this trajectory must 
be explicitly taken into account in the feedback controller 
design. Note that in literature pinning synchronizability 
is often called pinning controllability. Here we use the 
term synchronizability to avoid confusion with the clas¬ 
sical notion of controllability discussed in Secs. [IT] and 

EH 

A controlled network is described by 

N 

Xj = f(xi) - cr^gyh(xj) -P 6 ^Ui{t), (99) 

where <5^ = 1 for pinned nodes and 0 otherwise, and 


U,{t) = cr[p*(s(t)) - p,(Xi(f))] 


( 100 ) 


is the d-dimensional linear feedback controller ( Li et al.\ 
2004[ Wang and Chen 2002a), Pi(x(t)) is the pin¬ 
ning function that controls the input of node i, and 
s{f) is the desired synchronization trajectory satisfying 
s(t) = f(s(t)). Note that in the fully synchronized state 
Xi(t) = X 2 (t) = • • • = XAr(t) = s(t), we have Ui{t) = 0 
for all nodes. The form of the linear feedback controller 
(llOOl) implies that the completely synchronized state is a 


natural solution of the controlled network (99). 


Similar to spontaneous synchronization, we must de¬ 
rive the necessary and sufficient conditions for pinning 
synchronization. These conditions are more important 
from the control perspective, because they are the pre¬ 
requisite for the design of any practical controller. If we 
focus on the local (or global) stability of the synchro¬ 
nized manifold of the controlled network ( |99[ ), we ob¬ 
tain the necessary (or sufficient) condition for pinning 
synchronization, describing the local (or global) pinning 
synchronizability. 

Local pinning synchronizability: Given the presence of 
inhomogeneous dynamics at the controlled and uncon¬ 
trolled nodes, the MSF approach can not be directly 


applied to the controlled network (99). Instead, we 


first introduce a virtual node whose dynamics follows 
s{t) = f(s(t)), representing the desired synchronization 
solution ( Sorrentino et al.\ 2007b Zou and Chen[ 2008| ). 
The extended system now has A^-Pl nodes: yi{t) = Xi{t) 
for i = 1, • • • , A7; and yN+i{t) = s{t). The virtual node 
is connected to each pinned node. 

We choose the pinning function 


Pi(x) = K*h(x) 


( 101 ) 


with control gains Ki > 0, parameters that capture the 
relationship between the magnitude of h(x) and Pi(x). 
By defining the pinning function via (101) we can then 
rewrite (99) in the form of (90), with an effective coupling 
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matrix satisfying the zero row-sum condition, allowing us 


to apply the MSF approach. Indeed, plugging (1011 into 


(lOOl, we have Ui(t) = crH;i[h(s(t)) — h(xi(t))] and (99) 
becomes 


AT-1-1 

= f(y*) - CT ^ rn^jhiyj) 
1=1 


( 102 ) 


where 


M = 


9ii + SiKi 

912 

9lN 

-SlKi 

921 

922 + 62H2 

93N 

— 62 K2 

9ni 

9N2 

■ • • 9nn + Snkn 

—Snkn 

0 

0 

0 

0 


(103) 

is the effective coupling matrix of the (iV-l-l)-dimensional 
extended system. Apparently, M is a zero row-sum ma¬ 
trix, hence we can sort its eigenvalues as 0 = Ai < 
ReA 2 < ••• < RcAat+i. We can now apply the MSF 
approach to numerically explore the local stability of the 


synchronization manifold of the controlled network (1021. 


The role of the control gain (k^), coupling gain (cr), 
and the number and locations of the pinned nodes, on 
local pinning synchronizability has been systematically 
studied ( Sorrentino et al] |2007b| . Consider for ex¬ 
ample a Barabasi-Albert (BA) scale-free network of N 
identical Rossler oscillators coupled in x and z direc¬ 
tions. By assuming ki = ■ ■ ■ = kn = k, it was found 
that for a wide range of coupling gain a, the eigenra- 
tio = ReAAr-i-i/ReA 2 of the new coupling matrix 

M is minimized and hence the local pinning synchroniz¬ 
ability is maximized around a specific cr-dependent value 
of the control gain k. In other words, too large or too 
small control gain can reduce the network pinning syn¬ 
chronizability (Fig. [M^,b). In contrast, the number of 
pinned nodes, regardless if they are chosen randomly or 
selectively within the network, has a monotonic impact 
on pinning synchronizability: Controlling more nodes al¬ 
ways enhances the network pinning synchronizability, in 
line with our intuition (Fig. [Mj :,d). Furthermore, selec¬ 
tive pinning, when the nodes are chosen in the order of 
decreasing degree, yields better synchronizability than 
random pinning. 

Global pinning synchronizability. By describing the 


time evolution of the controlled network (102) in terms of 


the error dynamics, we can map the global pinning syn¬ 
chronizability of (1021 to the global asymptotic stability 


of the synchronized manifold, which can be studied via 
Lyapunov stability theory. 

If the desired asymptotic trajectory is an equilibrium 
point (s = f(s) = 0), we can derive sufficient conditions 
for globally stabilizing the pinning controlled network (Li 


et al.\ 2004). For a more general desired trajectory, it has 


been shown that a single feedback controller can pin a 



FIG. 36 (Color online) Local pinning synchronizability of 
scale-free networks. The local pinning synchronizability is 
quantified by the eigenratio = ReAiv+i/ReA 2 of the 

extended system (102l. The calculation was performed for 
N = 10® identical Rossler oscillators coupled in x and z 
directions, with coupling gain a and a p fraction of pinned 
nodes, placed on a Barabasi-Albert (BA) scale-free network 
with mean degree (fc) = 4. (a-b) We choose p = 0.1 fraction 
of nodes to pin and study the impact of control gain k on local 
pinning synchronizability with coupling gain a = 0.3 (a) and 
2.8 (b), respectively. We hnd that in both cases the eigen¬ 
ratio = ReAjv-i-i/ReA 2 of the new coupling matrix M 

is minimized and hence the local pinning synchronizability is 
maximized around a specific cr-dependent value of the control 
gain K. (c-d): We study the impact of the fraction of pinned 
nodes on local pinning synchronizability: (c) a = 0.3, k = 10. 
(d) a — 2.8, K = 1.5. The horizontal continuous lines (red) 
represent the eigenratio of the corresponding uncontrolled 
system (99). We find that the number of pinned nodes, re¬ 
gardless if they are chosen randomly or selectively within the 
network, has a monotonic impact on the pinning synchroniz¬ 
ability. Controlling more nodes always enhances the network 
pinning synchronizability. In all plots squares represent the 
case of random pinning, i.e., a p fraction of nodes is randomly 
chosen to be pinned. In (c) and (d), triangles represent the 
case of selective pinning, where nodes have been sorted in the 
order of decreasing degree and t he top p fraction of the n odes 
are chosen to be pinned. After (Sorrentino et al. 2007bI. 


complex network to a homogenous solution, without as¬ 
suming symmetry, irreducibility, or linearity of the cou¬ 
plings (Chen et al. 2007). 

If the oscillator dynamics f(x) fulfills 


f(zi) - f(z2) = J'zi,Z2(Zl - Z2), Vzi,Z2GK‘', ( 104 ) 

where is bounded, i.e., there exists a posi¬ 
tive constant a such that for any zi, Z 2 G < 

a, then we can derive tractable sufficient conditions for 
global pinning synchronizability in terms of the network 
topology, the oscillator dynamics, and the linear state 
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feedback (^ 

^orfiri and di Bernardo 

2008 ). Note that con- In other words, the control gain Ki varies in time and 

dition ( 

10 ^ 

:) applies to 

a large variety of chaotic oscilla- adapts to the error vector ei(t) = s(t) — Xi(t), that de- 

tors ( 

Jiang et al. 

2003 

1. The results indicate that for a scribes the deviation of the oscillator i from the reference 


connected network, even for a limited number of pinned 
nodes, global pinning synchronizability can be achieved 
by properly selecting the coupling strength and the feed¬ 


back gain (Chen et at. 2007). 


If h(x) = Fx and the oscillator dynamics f(x) satisfies 
(x - y)T [f (x, t) - f (y, t)] < (x - y)TKr(x - y) (105) 
for a constant matrix K, sufficient conditions for global 


pinning synchronizability can also be derived (Song and 
Cao 2010 lYu et al. 2013, 2009bI. Note that the con¬ 
dition (105) is so mild that many systems, from Lorenz 


system to Chen system, Lii system, recurrent neural net¬ 


works, Chua’s circuit satisfy this condition (Yu et al. 


2013). Counterintuitively, it was found that for undi¬ 


rected networks, the small-degree nodes, instead of hubs, 
should be pinned first when the coupling strength a is 
small ( |Yu et al. 2009b). For directed networks, nodes 
with very small in-degree or large out-degree should be 
pinned first ( Yu et aT\ |2013 ). This result can be under¬ 
stood by realizing that low in-degree nodes receive less 
information from other nodes and hence are less “influ¬ 
enced” by others. In the extreme case, nodes with zero 
in-degree will not be “influenced” by any other nodes, 
hence they must be pinned first. On the other hand, 
large out-degree nodes can influence many other nodes, 
hence it makes sense to pin them first. 


3. Adaptive pinning control 


Implementing the linear feedback pinning controller 


(1001 requires detailed knowledge of the global network 


topology. This is because we have to check whether there 
are possible coupling and control gains that ensure pin¬ 
ning synchronizability. Yet, in practice we do not always 
have access to the global network topology. Given this 
limitation, recently adaptive control has been proposed 
for pinning synchronization, in which case a controller 
adapts to a controlled system with parameters that vary 
in time, or are initially uncertain, without requiring a de¬ 


tailed knowledge of the global network topology (DeLellis 


eT^ 20TT1 |20T01 IWang 

et aL\ 200811 Wang and Slotinel 

2006 Wang et al. 

2010 

Zhou et al. 2008). As we dis- 


cuss next, many different strategies have been designed 
to tailor the control gains, coupling gains, or to rewire the 
network topology to ensure pinning synchronizability. 

(i) Adaptation of control gains: To adapt the control 


signal s(t). If the individual dynamics f(x) satisfies the 
Lipschitz condition, then the global stability of this adap¬ 
tive strategy can be assured. 

(ii) Adaptation of coupling gains: The coupling gain 
aij, defining the mutual coupling strength between node 


pair {i,j), can also be adapted using (jDeLellis et al. 


d-y(t) =?7y|ei(t)-ej(t)|^ 


(107) 


This strategy is very effective in controlling networks of 
quadratic dynamical systems, where the dynamics f(x, t) 
of each oscillator satisfies (x —y)'^[f(x, t) — f(y,t)] — (x — 
y)"^A(x — y) < —a;(x — y)'^(x — y). Here, A is an d x d 
diagonal matrix and w is a real positive scalar. 


Note that the adaptive strategies (106) and (107) are 


based on the local error vectors of nodes or between 
neighboring nodes, hence they avoid the need for a prior 
tuning of the control or coupling gains. This is attractive 
in many circumstances. However, these adaptive strate¬ 
gies still require a prior selection of the pinned nodes 
based on some knowledge of the network topology. This 
limitation can be avoided by choosing pinned nodes in 
an adaptive fashion, as we discuss next. 

(hi) Adaptive selection of pinning nodes: Adaptive pin¬ 
ning can be achieved by assuming the pinning node in¬ 
dicator 5i to be neither fixed nor binary. A common 
approach is to introduce 

5i{t) = bi{t) (108) 

where hi{t) satisfies the dynamics 

H-TT—^=5(1®*!)- (109) 

Abi 

In other words, bi{t) follows the dynamics of a unitary 
mass in a potential U{bi) subject to an external force g 
that is a function of the pinning error and a linear 
damping term described by C,bi. This is termed as the 
edge-snapping mechanism. For convenience, U{-) can be 
chosen as a double-well potential: U[z) = kz^{z — 1)^, 
where the parameter k defines the height of the barrier 


between the two wells. Then (109) has only two stable 


equilibria, 0 and 1, describing whether node i is pinned or 
not, respectively. Sufficient conditions for the edge snap¬ 


ping mechanism (109) to drive the network to a steady- 


state pinning configuration have been derived (DeLellis 


gain Ki (1011, representing the ratio between the pinning ^OH ). The key advantage of the adaptive selection 


function and output function, we choose the control in¬ 
put Ui(f) = —6iKi{t){xi{t) — s), and the control gains 
as (Wang et al. 2008[ Zhou et at. 2008) 


ki{t) = qi\ei{t)\. 


(106) 


of pinning nodes is that we don’t have to choose the nodes 
we need to pin before we design the controller. Instead, 
we can select them as we go in an adaptive fashion. 

(iv) Adaptation of the network topology: We can ensure 
synchronization by adapting the network topology. Spe- 
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cially, we can set each off-diagonal element of the Lapla- 
cian matrix of the network as 




( 110 ) 


and its variants can be interpreted as decentralized feed¬ 
back control system with time-varying network structure, 
offering a better understanding of the origin of collective 


where is the mutual coupling strength between 


behavior ( 

Jadbabaie et al. 

2003 

Moreau 

2005 

Olfati- 

Saber 

2006 

Ren and Beard 

2005 

I. 


node pair (i, j), which is adapted as in (107). The weight 


aij(t) is associated to every undirected edge of the target i. Vicsek Model and the Alignment Problem 
pinning edge and is adapted as 


The Vicsek model explains the origin of alignment, a 


da,: 


= c(|ey|), i,j = 1,... j, 

( 111 ) 

where eij{t) = ej{t) — and U{-) can be again cho¬ 

sen as a double-well potential so that ( [ 111 ) has only two 
stable equilibria, 0 and 1. In this case, the target net¬ 
work topology evolves in a decentralized way. The local 
mismatch of the trajectories can be considered as an ex¬ 


key feature of flocking behavior (Vicsek et al. 1995). It 


ternal forcing on the edge dynamics ( 111 ), inducing the 


activation of the corresponding link, i.e. Uij = 1 . 

The above adaptive strategies cope better when pin¬ 
ning controllability using a non-adaptive or static ap¬ 
proach is initially not feasible. They are also successful in 
ensuring network synchronization in the presence of per¬ 


turbations or deterioration, like link failures (Jin et al. 


2012 ). 


Taken together, we have multiple strategies to force a 
networked system to synchronize. The discussed tools 
have a wide range of applications for systems in which 
a synchronized state is desired. In some cases synchro¬ 
nization can be harmful, like in the case of synchronized 
clients or routers that cause congestion in data traffic on 


the Internet (Li and Chen, 2003), or in schizophrenia. In 


this case the synchronized state can be destroyed by the 


addition of a single link with inhibitory coupling (Slotine 


et al. 2004). 


B. Flocking of multi-agent dynamic systems 

The flocking of birds, shoaling of fish, swarming of in¬ 
sects, and herding of land animals are spectacular man¬ 
ifestations of coordinated collective behavior of multi¬ 
agent systems. These phenomena have fascinated sci¬ 
entists from diverse disciplines, from ecologists to physi¬ 


cists, social and computer scientists (Olfati-Saber 2006 


Vicsek and Zafeiris 2012). Many models have been pro¬ 


posed to reproduce the behavior of such self-organized 
systems. The first widely-known flocking simulation was 
primarily motivated by the visual appearance of a few 
dozen coherently flying objects, e.g., imaginary birds and 


spaceships (Reynolds 1987). Yet, the quantitative inter¬ 


pretation of the emerging behavior of huge flocks in the 
presence of perturbations was possible only following the 
development of a statistical physics-based interpretation 


of flocking obtained through the Vicsek model (Vicsek 


is a discrete-time stochastic model, in which autonomous 
agents move in a plane with a constant speed vq, initially 
following randomly chosen directions. The position of 
agent i changes as 


Xi(t -k 1) = X,(t) -k v,(t -k 1), 


( 112 ) 


where the velocity of each agent has the same absolute 
value vq. The direction of agent i is updated using a local 
rule that depends on the average of its own direction and 
the directions of its “neighbors”, i.e. all agents within a 
distance r from agent i (FigjS^. In other words. 


9i{t -k 1) — {0i(t))r -k Ai(t). 


(113) 


Here {9i[t))r = arctan [(sin0(t))r/(cos0(t))r] denotes the 
average direction of the agents (including agent i) within 
a circle of radius r. The interaction radius r can be set as 
the unit distance, r = 1. The origin of the alignment rule 


(1131 can be the stickiness of the agents, hydrodynamics, 
could be pre-programmed, or based on information pro¬ 


cessing (Vicsek and Zafeiris 2012). The perturbations 


are contained in Ai(t), which is a random number taken 
from a uniform distribution in the interval [— 77 / 2 , 77 / 2 ]. 
Therefore the final direction of agent i is obtained after 
rotating the average direction of the neighbors with a ran¬ 
dom angle. These random perturbations can be rooted 
in any stochastic or deterministic factors that affect the 
motion of the flocking agents. 

The Vicsek model has three parameters: (i) the agent 
density p (number of agents in the area L^); (ii) the 
speed Vq and (iii) the magnitude of perturbations 77 . The 
model’s order parameter is the normalized average veloc¬ 
ity 


1 


Nvq 


N 

E' 

i=l 


(114) 


For small speed WQ) if we decrease the magnitude of 
perturbations 77 , the Vicsek model displays a continuous 
phase transition from a disordered phase (zero average 
velocity (f), implying that all agents move independently 
of each other. Fig. [3^) to an ordered phase when almost 
all agents move in the same direction, through a spon¬ 
taneous symmetry breaking of the rotational symmetry 


et al. 1995). As we discussed next, the Vicsek model (Fig. 381). This much studied kinetic phase transition 
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FIG. 37 (Color online) Vicsek model. The direction of agent i 
at time t + 1 (shown in red) is the average of its own direction 
and the directions of all other agents at a distance less than 
r to agent i at time t (shown in grey). Agents outside this 
circle (shown in white), do not contribute to the direction of 
agent i at time t+1. 


takes place despite the fact that each agent’s set of near¬ 
est neighbors change with time as the system evolves and 
the absence of centralized coordination. 

Numerical results indicate that the phase transition is 
second-order and the normalized average velocity cj) scales 
as 


^ ^ [Vc{p) - vf, 


(115) 


where the critical exponent /? « 0.45 and r]c{p) is the crit¬ 
ical noise for L —oo (Vicsek et a/., [1995 ). Many studies 
have explored the nature of the above phase transition 
(whether it is first or second order), finding that two fac¬ 
tors play an important role: (i) the precise way that the 
noise is introduced into the system; and (ii) the speed vq 
with which the agents move ( [Aldana et al. [2007 
Baglietto and Albano[ |2009[ [Gregoire and Chate 
Pimentel et al. 20081. 


2009 


2004 


The Vicsek model raises a fundamental control prob¬ 
lem: Under what conditions can the multi-agent system 
display a particular collective behavior? Behind each 
flock of collectively moving agents, like biological organ¬ 
isms or robots, there is a dynamically changing or tem¬ 
poral network, where two agents are connected if they 
interact, e.g. if their distance is under a certain thresh¬ 
old. Since the agents are moving, the network of momen¬ 
tarily interacting units evolves in time in a complicated 
fashion. 

To offer a control theoretical explanation for the emer¬ 
gence of the ordered phase in the Vicsek model, we con¬ 


sider the following updating rule (Jadbabaie et al. 20031: 


6i{t -|-1) — 


1 


1 -|- ki{t) 


0^{t) 


E 


(116) 


Though the scalar average in (116) is fundamentally dit 


ferent from the vectorial average in (113), this updating 


rule still captures the essence of the Vicsek model in the 



FIG. 38 Emergence of order in the Vicsek Model. The panels 
show the agent velocity for varying values of the density and 
the noise level. The actual velocity of an agent is indicated by 
a small arrow while their trajectories for the last 20 time steps 
are shown as short continuous curves.The number of agents 
is A = 300, and the absolute velocity is vq = 0.03. (a) At t = 
0, the positions and the direction of velocities are randomly 
distributed. L = 7, »7 = 2.0 (b) For small densities (L = 25) 
and noise {rj — 0.1) level, the agents form groups that move 
together in random directions, (c) At higher densities {L = 7) 
and noise {rj = 2.0) the agents move randomly with some 
correlation, (d) When the density is large (L = 5) and noise is 
small (t; = 0.1), the motion becomes ordered on a macroscopic 
scale and all agents tend to move in the same spontaneously 
selected direction. After (Vicsek et al. 19951. 


absence of perturbation. More importantly, (116) can be 
considered as a decentralized feedback control system 


e{t + i) = e{t) + u{t) (117) 

with the control input 

u(t) = -(D,p)+I)-iL,p)e(t). (118) 

Here Lp = Dp — Ap is the Laplacian matrix of graph Gp 
with p € V. A.p is the adjacency matrix of graph Gp and 
Dp is a diagonal matrix whose ith diagonal element is 
the degree of node i in the graph Gp. a{t) : 0,1, •••—>■ 7^ 
is a switching signal whose value at time t is the index of 
the interaction graph at time t, i.e., G{t). 

If r is small, some agents/nodes are always isolated, 
implying that Git) is never connected. If r is large, then 
G{t) is always a complete graph. The situation of inter¬ 
est is between the two extremes. The goal is to show that 
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despite the fact that each agent’s set of nearest neigh¬ 
bors changes in time. These control theoretical results 
suggest that to understand the effect of additive noise, 
we should focus on how noise inputs effect connectivity 
of the associated neighbor graphs. For example, the nu¬ 
merical finding that, for a fixed noise beyond a critical 
agent density all agents eventually become aligned, can 
be adequately explained by percolation theory of random 
graphs (Jadbabaie et al.. 20031. 


FIG. 39 Kinetic phase transition in the Vicsek model, (a) 
The normalized average velocity (</>) versus the magnitude of 
perturbations (noise rj) in cells of various sizes (L) with a 
fixed density p = N/L'^ = 0.4. As rj decreases, (f> increases, 
implying the emergence of order in the Vicsek model, (b) 
Dependence of 4> on [pc(T) — pj/pc in log-log scale. The slope 
of the lines is associated with the critical exponent j3 for which 
we get j3 = 0.45 ± 0.07. The scaling behavior of 0 observed 
in such a kinetic phase transition is analogous to what we 
often observe in continuous phase transitions in equilibrium 


systems. After (Vicsek et al. 19951. 


for any initial set of agent directions 0(0) and for a large 
class of switching signals the directions of all agents will 
converge to the same steady state 0ss, reaching align¬ 
ment asymptotically. Mathematically, this means that 
the state vector 9{t) converges to a vector of the form 
0ssl with 0SS the steady state direction, i.e., 


lim 6{t) = 9s, 

t—^OC) 


(119) 


where 1 = (1, • • • , l)^xi’ representing the case when all 
agents move in the same direction. 

If G{t) is connected for all t > 0, then we can prove 


that alignment will be asymptotically reached (Jadbabaie 


et al.. 20031. But this condition is very stringent. It can 


be relaxed by considering that the agents are linked to¬ 
gether across a time interval, i.e., the collection or union 
of graphs encountered along the interval is connected. It 
has been proven that if the N agents are linked together 
for each time interval, then the alignment will be asymp¬ 
totically reached (Jadbabaie et al. [2003 ). This result has 


been further extended by proving that if the collection of 
graphs is ultimately connected, i.e., there exists an ini¬ 
tial time to such that over the infinite interval [tg, oo) the 


ment is asymptotically reached ( 

Moreau 

2005). 

Though the control theoretical analysis ( 

Jadbabaie 

et al. 

2003 

Moreau 

2005 

Ren and Beard 2005 

) is deter- 


ministic, ignoring the presence of noise, it offers rigorous 
theoretical explanations, based on the connectedness of 
the underlying graph, for some fundamental aspects of 
the Vicsek model. For example, by applying the nearest 
neighbor rule, all agents tend to align the same direc¬ 
tion despite the absence of centralized coordination and 


2. Alignment via pinning 


While the virtue of the Vicsek model is its ability to 
spontaneously reach an ordered phase, we can also ask 
if such a phase can be induced externally. Therefore, we 
consider an effective pinning control strategy in which 
a single pinned node (agent) facilitates the alignment of 
the whole group. This is achieved by adding to the Vic¬ 
sek model an additional agent, labeled 0, which acts as 
the group’s leader. Agent 0 moves at the same constant 
speed Vo as its N followers but with a fixed direction 0o, 
representing the desired direction for the whole system. 
Each follower’s neighbor set includes the leader whenever 
it is within the follower’s circle of radius r. Hence we have 


9i{t+l) — 


1 + ki{t) + bilt) 


9j{t) + bi(t)9o 

( 120 ) 

where bi(t) = 1 whenever the leader is a neighbor of 
agent i and 0 otherwise. It has been proved that if the 
{N + 1) agents are linked together for each time inter¬ 


val, then alignment will be asymptotically reached (Jad- 


babaie et al.\ 20031. In other words, if the union of graphs 


of the {N -I-1) agents encountered along each time inter¬ 
val is connected, then eventually all the follower agents 
will align with the leader. 


3. Distributed flocking protocols 


Alignment, addressed by the Vicsek model, is only one 
component of flocking behavior. Indeed, there are three 
heuristic rules for flocking (Reynolds 1987): (i) Cohe¬ 


sion: attempt to stay close to nearby flockmates; (ii) 
Separation: avoid collisions with nearby flockmates; and 
(iii) Alignment: attempt to match velocity with nearby 
flockmates. 

We therefore need a general theoretical framework to 
design and analyze distributed flocking algorithms or pro¬ 
tocols that embody these three rules. The formal ap¬ 
proach described next extracts the interaction rules that 


can ensure the emergence of flocking behavior (Olfati- 


Saber 2006). 


Consider a gradient-based flocking protocol equipped 
with a velocity consensus mechanism, where each agent 
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FIG. 40 Geometry of flocking and fragmentation, (a) Lattice- 
type flocking configuration in D = 2. In this ideal case, each 
agent is at the same distance from all of its neighbors on the 
proximity graph, (b) A quasi-lattice for D = 2 with N = 150 
nodes, (c) Fragmentation phenomenon, where agents merge 
form a few groups and different groups are moving in differ¬ 
ent directions. This configuration will never lead to flocking 
behavior. After (Olfati-Saber 20061. 


is steered by the control input 

u, =ff + ff. (121) 

The first term 

ff = (122) 

is gradient-based and regulates the distance between 

agent i and its neighbors, avoiding the collision and co¬ 
hesion of the agents. This term is derived from a smooth 
collective potential function 14 (q), which has a unique 
minimum when each agent is at the same distance from 
all of its neighbors on the proximity graph G(q), repre¬ 
senting the ideal case for flocking. The second term 

ff = ^ a„(t)(pj-Pi) (123) 

leAfdb 


regulates the velocity of agent i to match the average 
velocity of its neighbors, being responsible for the ve¬ 
locity alignment. Here the weighted spatial adjacency 
matrix A(t) = [aij(t)] is calculated from the proximity 
network G(q). The flocking protocol (121 1 embodies all 
three rules of Reynolds. However, for a generic initial 
state and a large number of agents (e.g., N > 100), the 


protocol (1211 leads to fragmentation, rather than flock¬ 


ing (Olfati-Saber 2006), meaning that the agents sponta¬ 


neously form several groups, where different groups move 


in different directions (Fig. 40:). To resolve this fragmen¬ 


tation issue, we introduce a navigational feedback term 
to the control input of each agent 


u, = ff + ff 


(124) 


where 


f7 =-ci(qi - q 7 ) - C 2 (pi - p^) (125) 

drives agent i to follow a group objective. The group 
objective can be considered as a virtual leader with the 
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FIG. 41 (Golor online) Flocking behavior in mutli-agent sys¬ 
tems. After the application of the flocking algorithm (1241 
for a few seconds, the flocking of A = 100 agents in 2D is 
observed. After (Olfati-Saber 20061. 


following equation of motion: 


qj = p.y 

^P7 =f7(q7:P7)’ 

pD 


(126) 


where q 7 ,P 7 ,f 7 (q 7 ,P 7 ) S are the position, velocity, 
and acceleration (control input) of the virtual leader, re¬ 
spectively. By taking into account the navigational feed¬ 
back, the protocol ( |124 ) enables a group of agents to 
track a virtual leader that moves at a constant velocity. 


and hence leads to flocking behavior (Olfati-Saber 2006 1 . 


Note that protocol (124) requires all agents to be in¬ 


formed, i.e., to know the group objective, or equivalently, 
the current state (q 7 ,P 7 ) of the virtual leader. It turns 
out that this is not necessary for flocking. Motivated by 
the idea of pinning control, it has been shown that, even 
when only a fraction of agents are informed (or pinned), 
the flocking protocol (1241 still enables all the informed 
agents to move with the desired constant velocity. An 
uninformed agent will also move with the desired veloc¬ 
ity if it can be influenced by the informed agents from 
time to time (Su et at. 2009). Numerical simulations 


suggest that the larger the informed group is, the bigger 
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fraction of agents will move with the desired velocity (Su 

eFZj[200^ . 

If the virtual leader travels with a varying velocity 


p-y(t), the flocking protocol (124) enables all agents to 


eventually achieve a common velocity. Yet, this common 
velocity is not guaranteed to match p.^ (t). To resolve this 
issue, we can incorporate the acceleration of the virtual 


leader into the navigational feedback (125) as follows 


f? = f7(q7>P7) - Ci(qi - q.^) - C2(Pi - p.y). (127) 

The resulting protocol enables the asymptotic tracking 
of the virtual leader with a varying velocity, ensuring 
that the position and velocity of the center of mass of all 
agents will converge exponentially to those of the virtual 
leader (ISu et al.\ 120091). 


In summary, the combination of control theoretical and 
network science approaches can help us understand the 
emergence of order in multi-agent systems. These tools 
are indispensable if we wish to understand how to induce 
order externally, aiming to control the collective behavior 
of the system. 


VII. OUTLOOK 

Given the rapid advances in the control of complex 
networks, we have chosen to focus on a group of results 
that will likely stay with us for many years to come. The 
process of organizing the material has also exposed obvi¬ 
ous gaps in our knowledge. Therefore, next we highlight 
several research topics that must be addressed to realize 
the potential of the control of complex systems. Some 
of these may be addressed shortly, others, however, may 
continue to challenge the community for many years to 
come. 


A. Stability of Complex Systems 


Stability is a fundamental issue in the analysis and the 
design of a control system, because an unstable system 
is extremely difficult and costly to control, and such a 


system can also be potentially dangerous (Chen 2001 


Slotine and Li 1991). Loosely speaking, a system is sta¬ 


ble if its trajectories do not change too much under small 
perturbations. 

The stability of a nonlinear dynamical systems x = 
f (x, t) can be analyzed by the Lyapunov Stability The¬ 
ory (LST), without explicitly integrating the differential 
equation. LST includes two methods: (i) The indirect (or 
linearization) method, concerned with small perturbation 
around a system’s equilibrium points x* and the stability 
conclusion is inferred from a linear approximation of the 
nonlinear systems around this equilibrium point. This 
justifies the use of linear control for the design and analy¬ 
sis of weakly nonlinear systems, (ii) The direct method is 


based on the so-called Lyapunov function— an “energy¬ 
like” scalar function whose time variation can be viewed 
as “energy dissipation”. It is not restricted to small per¬ 
turbations and in principle can be applied to any dy¬ 
namical system. Yet, we lack a general theory to find a 
suitable Lyapunov function for an arbitrary system. We 
have to rely on our experience and intuition to formulate 
Lyapunov functions, like exploiting physical properties 


(such as energy conservation) and physical insights (Slo¬ 


tine and Li, 1991). 


For a wide range of complex systems certain diagonal- 
type Lyapunov functions are useful for stability analy¬ 
sis (Kaszkurewicz et al. 2000). More importantly, in 
many cases the necessary and sufficient conditions for the 
stability of nonlinear systems are also the necessary and 
sufficient conditions for the diagonal stability of a certain 
matrix associated to the nonlinear system. This matrix 
naturally captures the underlying network structure of 
the nonlinear dynamical system. 

Matrix diagonal stability is a well-known notion in sta¬ 
bility analysis since its introduction by Volterra around 
1930 in the context of ecological systems (Volterra 1931). 
Yet, its usefulness is limited by the difficulty of char¬ 
acterizing the class of large diagonally stable matrices. 
Though there are efficient optimization-based algorithms 
to numerically check if a given matrix is diagonally sta¬ 
ble (Boyd et al. 1994), there are no effective theoreti¬ 
cal tools to characterize general large diagonally stable 
matrices. Recently, however, necessary and sufficient di¬ 
agonal stability conditions for matrices associated with 
special interconnection or network structures were stud¬ 
ied (Arcak 2011 Arcak and Sontag 2008,2006), improv¬ 
ing our understanding of the stability of gene regulatory 
and ecological networks. More research is required to un¬ 
derstand stability, an important prerequisite for control. 

The stability concepts we discussed above consider per¬ 
turbations of initial conditions for a fixed dynamical sys¬ 
tem. There is another important notion of stability, i.e. 
structural stability, which concerns whether the qualita¬ 
tive behavior of the system trajectories will be affected by 


small perturbations of the system model itself (Andronov 


and Pontryagin 

1937 

Kuznetsov 

o 

o 


To formally define structural stability, we introduce the 
concept of topologically equivalence of dynamical sys¬ 
tems. Two dynamical systems are called topologically 
equivalent if there is a homeomorphism h : 
mapping their phase portraits, preserving the direction 
of time. Consider two smooth continuous-time dynamical 
systems (1) x = f(x); and (2) x = g(x). Both (1) and (2) 
are defined in a closed region D G (see Fig. 42). Sys¬ 
tem (1) is called structurally stable in a region Dq C D ii 
for any system (2) that is sufficiently C^-close to system 
(1) there are regions U,V C D, and Dq C U, Dq C V 
such that system (1) is topologically equivalent in U 
to system (2) in V (see Fig. 42 1 ). Here, the systems 
(1) and (2) are C^-close if their “distance”, defined as 
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FIG. 42 (Color online) Structural stability. (a) An¬ 
dronov’s definition of structural stability, (b-d) Phase por¬ 
traits of structurally unstable planar systems. This figure is 
redrawn from Figures 2.19 and 2.20 of ( Kuznetsov[ 20041. 


df(x) 

dx 


dg(x 


dx 


}■ 


is small 


di = sup^g£, |||f(x) -g(x)|| -h 
enough. 

For a two-dimensional continuous-time dynamical sys¬ 
tem, the Andronov-Pontryagin criterion offers sufficient 


and necessary conditions for structural stablility (An 
dronov and Pontryagin[ 19371. A smooth dynamical sys- 
tern X = f(x),x S is structurally stable in a region 
Dq C if and only if (i) it has a finite number of equi¬ 
librium points and limit cycles in Dq, and all of them are 
hyperbolic; (ii) there are no saddle separatrices return¬ 


ing to the same saddle (see Fig. 42 b,c) or connecting two 
different saddles in Dq (see Fig. |42[1). It has been proven 
that a typical or generic two-dimensional system always 
satisfies the Andronov-Pontryagin criterion and hence is 


structurally stable (Peixoto 1962). In other words, struc¬ 


tural stability is a generic property for planar systems. 
Yet, this is not true for high-dimensional systems. 

For A^-dimensional dynamical systems, Morse and 
Smale established the sufficient conditions of structural 


stability (Smale 1961 1967). Such systems, often called 


Morse-Smale systems, have only a finite number of equi¬ 
librium points and limit cycles, all of which are hyper¬ 
bolic and satisfy a transversaility condition on their sta¬ 
ble and unstable invariant manifolds. 

The notion of structural stability has not been well 
explored in complex networked systems. 


B. Controlling Adaptive Networks 


namic networks in control theory (Mesbahi 2005 Mes- 


bahi and Egerstedt, 2010), are collections of units that 


interact through a network, whose topology evolves as 
the state of the units changes with time. Adaptive net¬ 
works are a special class of temporal networks, whose 


edges are not continuously active ( 

Holme and Saramaki 

2012 

Karsai et aLl |2011[ |Pan anc 

Li 

2014 

Posfai and 


Hovel 


2014). If the temporal order of the network snap¬ 


shots at different time points depend on the states of 
the nodes, then the temporal network is adaptive. A 
special case of adaptive networks are switched systems, 
which consist of a family of subsystems and a switching 


law that orchestrates the switching among them (Xie and 


Wang 2003 Xie et al. 2002). For switching systems, we 


can design the switching signal among different subsys¬ 
tems and hence the switching law may be independent 
from the states of the nodes. 

Mycelial fungi and acellular slime molds grow as self- 
organized networks that explore new territory for food 
sources, whilst maintaining an effective internal trans¬ 
port system to resist continuous attacks or random dam¬ 


age (Fessel et ai. 2012). Honed by evolution, these bi¬ 


ological networks are examples of adaptive transporta¬ 
tion networks, balancing real-world compromises be¬ 


tween search strategy and transport efficiency (Tero 


et ai. 2010). 


The genome is also an intriguing example of an adap¬ 
tive network, where the chromosomal geometry directly 
relates to the genomic activity, which in turn strongly 
correlates with geometry ( [Rajapakse et al. . 2011). Simi¬ 
larly, neuronal connections (synapses) in our brains can 
strengthen or weaken, and form in response to changes 
in brain activity, a phenomenon called synaptic plastic¬ 


ity (Bayati and Valizadeh, 2012 Perin et al. 2011). 


A comprehensive analytical framework is needed to 
address the control of adaptive, temporal and co¬ 
evolutionary networks. This framework must recognize 
the network structure itself as a dynamical system, to¬ 
gether with the nodal or edge dynamics on the network, 
capturing the feedback mechanisms linking the structure 
and dynamics. Studying the controllability of such sys¬ 
tems would be a natural starting point because seemingly 
mild limitations on either the network structure or the 
dynamical rules may place severe constraints on the con¬ 


trollability of the whole system (Rajapakse et al. 2011). 


Identifying these constraints is crucial if we want to re¬ 
frain from improving systems that already operate close 
to their fundamental limits. 


Adaptability, representing a system’s ability to re¬ 
spond to changes in the external conditions, is a key 
characteristic of complex systems. Indeed, the structure 
of many real networks co-evolves with the dynamics that 
takes place on them, naturally adapting to shifting envi¬ 


ronments (Gross and Sayama 2009). 


C. Controlling Networks of Networks 

Many natural and engineered systems are composed 
of a set of coupled layers or a network of subsystems, 
characterized by different time scales and structural pat- 


Adaptive networks, also known as state-dependent dy- terns. New notions, from multiplex networks (Boccaletti 
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et al. 

2014 

Kivel et al. 2014 

) to networks of net- affected by large fluctuations due to the low copy number 

works ( 

D’Agostino and Scala 

2014 

Gao et al. 2014a 

1, of transcription factors. The extrinsic noise of a biolog- 


have been recently proposed to explore the properties 
of these systems, focusing mainly on their structural in¬ 
tegrity and robustness. Consider a multiplex network, 
i.e. a set of coupled layered networks, whose different 
layers have different characteristics. We can model such 
a system as a layered network, whose interconnections 
between layers capture the interactions between a node 
in one layer and its counterpart in another layer. Simi¬ 
larly, in a network of networks each node itself is a net¬ 
work or a multi-input/multi-output (MIMO) subsystem. 
Different nodes/subsystems could have totally different 
dimensions and dynamics. This is rather different from 
the control framework discussed in much of this paper, 
where we typically assumed that all the nodes share the 
same type of dynamics or even just scalar dynamics (with 
state variables xt gR for all nodes). 

Developing a framework to control networks of net¬ 
works is a necessary step if we wish to understand the 
control principles of complex systems. Early attempts 
have focused on the issues of controllability or observ- 


ability with linear dynamics (Chapman et al. 2014 

Menichetti et al.. |2015[ 

Wang et al. 

2015 Yuan et al. 

2014 Zhang et al. 

2016 

Zhou 2015). 

For example, some 


controllability conditions on the overall network topol¬ 
ogy, the node dynamics, the external control inputs and 
the inner interactions have been derived for a networked 
MIMO system (Wang et al. 2015). Interestingly, the 


controllability of the networked MIMO system is an in¬ 
tegrated result of multiple factors, which cannot be de¬ 
coupled into the controllability of the individual subsys¬ 
tem or the properties solely determined by the network 
topology. Despite these efforts, we lack a general frame¬ 
work to systematically explore the control of networks of 
networks. Yet, the problem’s importance will likely trig¬ 
ger more research in both network science and control 
theory. 


D. Noise 

Complex systems, especially biological systems, are 
noisy. They are affected by two kinds of noise: the in¬ 
trinsic randomness of individual events and the extrinsic 
influence of changing environments ([Hilfinger and Pauls- 


son 


2011 


Lestas et a/., 2010). Consider, for example, 
regulatory processes in a cell. The intrinsic noise is 
rooted in the low copy number of biomolecules or diffu¬ 
sive cellular dynamics. In particular, if N is the number 
of molecules in the system, fluctuations in N lead to sta¬ 
tistical noise with intensity in the order of For 

large N, we can assume that a continuous deterministic 
dynamics effectively describes the changes of the average 
concentrations. However, for small N the statistical noise 
cannot be ignored. For example, gene regulation may be 


ical system is mainly due to the changing environments 
experienced by the system. The environmental change 
may have microscopic origin (like cellular age/cell cycle 
stage and organelle distributions) or can be related to 
the macroscopic physical or chemical environment (like 
illumination conditions, temperature, pressure and pH 
level). To infer or reconstruct the states of a biological 
system, we also need to deal with the measurement er¬ 
ror, which is independent of the biological system and 
can also be considered as extrinsic noise. 

Both internal and external noises are known to affect 
the control of complex systems. At this time we lack 
a full understanding on the role of noise or stochastic 
fluctuations on the control of complex systems. 


E. Controlling Quantum Networks 

Quantum control theory aims to offer practical meth¬ 
ods to control quantum systems. Despite recent progress, 


quantum control theory is still in its infancy (Dong and 


Petersen, 2010), for several reasons. First, in classical 


control it is assumed that the measurement does not af¬ 
fect the measured system. In contrast, in quantum con¬ 
trol it is difficult, if not impossible, to acquire information 
about quantum states without destroying them. Second, 
some classes of quantum control tasks, like controlling 
quantum entanglement and protecting quantum coher¬ 
ence, are unique for quantum systems. In other words, 
there are no corresponding tasks in classical control the¬ 
ory. 

The notion of quantum networks has been recently 


proposed by the quantum information community (Acm 

et a/.l 120071 

Cuquet and Calsamig 

ia 2009 2012 

Czekaj 

et al.. 2012 

Hahn et al. 2012 

Lapeyre et al. 

2009 

Perseguers 

1010 Perseguers et al. 

2008,2013 2010), of- 


fering fresh perspectives in the field of complex networks. 
In a quantum network, each node possesses exactly one 
qubit for each of its neighbors. Since nodes can act on 
these qubits, they are often called “stations”. The edge 
between two nodes represents the entanglement between 
two qubits. The degree of entanglement between two 
nodes can be considered as the connection probability (p) 
in the context of classical random graphs. 

In a classical random graph if we let p scale with 
the graph size as p ~ increasingly complex sub¬ 
graphs appear as z exceeds a series of thresholds. For 
example, for z < —2 almost all graphs contain only iso¬ 
lated nodes and edges. When z passes through —3/2 (or 
—4/3), trees of order 3 (or 4) suddenly appear. As z 
approaches —1, trees and cycles of all orders appear (Al¬ 
bert and Barabasi [2002 ). Surprisingly, in quantum net¬ 
works any subgraph can be generated by local operations 
and classical communication, provided that the entangle- 
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merit between pairs of nodes scales with the graph size 
as p ~ N~'^ (Perseguers et al. 2010| . In other words, 
thanks to the superposition principle and the ability to 
coherently manipulate the qubits at the stations, even for 
the lowest non-trivial connection probability that is just 
sufficient to get simple connections in a classical graph, 
we obtain quantum subgraphs of any complexity. 

This result illustrates that quantum networks have 
unique properties that are impossible in their classical 
counterparts. Hence, the control of quantum complex 
networks will require new methodologies. 


F. Conclusion 

Revealing the control principles of complex networks 
remains a challenging problem that, given its depth and 
applications, will probably engage multiple research com¬ 
munities for the next decade. In this review we aimed 
to summarize in a coherent fashion the current body of 
knowledge on this fascinating topic. This forced us to 
explore key notions in control theory, like controllability 
and observability, but also to explore how to steer a com¬ 
plex networked system to a desired final state/trajectory 
or a desired collective behavior. There are many out¬ 
standing open questions to be addressed, advances on 
which will require interdisciplinary collaborations. We 
hope that this review will catalyze new interdisciplinary 
approaches, moving our understanding of control forward 
and enhancing our ability to control complex systems. 
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